i/2 


flD-R134  692  LASER  DIAGNOSTIC  DEVELOPMENT  AND  MEASUREMENT  AND 
'  MODELING  OF  TURBULENT  FL  (U)  DAVTON  UNIV  OH  RESEARCH 

INST  L  KRISHNAHURTHV  ET  AL.  JUN  82  UDR-TR-82-04 
UNCLASSIFIED  AFWAL-TR-83-2044-PT-2  F33615-78-C-2005  F/G  20/4  NL 


:« 


r« 


m 


>- 

Cu 

o 


.  u 


IT  ■? 

r~: 


AFWAL-TR- 3 3-2044 
PART  II 


-/)  /2V  69  X 


LASER  DIAGNOSTIC  DEVELOPMENT  AND  MEASUREMENT  AND 
MODELING  OF  TURBULENT  FLOWF I  ELDS  OF  JETS  AND  WAKES 

Part  II 

Numerical  Predictions  of  Isotnermal  Flowfields 
in  a  Ducted  Centerbody  Combustor 


L.  Kr ishnamurthy ,  S.  0.  Park,  D.  J.  Wahrer,  and  H.  S.  Coch 
University  of  Dayton 
Research  Institute 
Davton,  Ohio  45469 


June  1983 


FINAL  REPORT  FOR  PERIOD  1  APRIL  1978  -  30  SEPTEMBER  1982 


Approved  for  Public  Release,  Distribution  Unlimited 


AERO  PROPULSION  LABORATORY 

AIR  FORCE  WRIGHT  AERONAUTICAL  LABORATORIES 
AIR  FORCE  SYSTEMS  COMMAND 

WRIGHT-PATTERSON  AIR  FORCE  BASE,  OH  45433 


NOT ICC 


Nhen  Government  drawings,  specifications,  or  other  data  are 
used  for  any  purpose  otner  than  in  connection  with  a  definitely 
related  Government  procurement  operation,  the  united  States 
Government  thereby  incurs  no  responsibility  nor  any  ooligation 
whatsoever;  and  the  fact  that  the  government  may  have  formulated, 
furnished,  or  in  any  way  supplied  the  said  drawings,  specifica¬ 
tions,  or  other  data,  is  not  to  be  regarded  by  implication  or 
otnerwise  as  in  any  manner  licensing  the  holder  or  any  other  per¬ 
son  or  corporation,  or  conveying  any  rights  or  permission  to 
manufacture  use,  or  sell  any  patented  invention  that  may  in  any 
way  be  related  thereto. 

This  report  has  been  reviewed  by  the  Office  of  Public  Affairs 
(ASD/PA)  and  is  releasable  to  the  National  Technical  Information 
Service  (NTIS).  At  NTIS,  it  will  be  available  to  the  general 
public,  including  foreign  nations. 

This  technical  report  has  been  reviewed  and  is  approved  for 
publication . 


u).  m.  _ 

W .  M .  ROQUEM^RE 
Fuels  Brancn 

Fuels  and  Lubrication  Division 
Aero  Propulsion  Laboratory 


FOR  THE  COMMANDER 


BENITO  P.  BOTTERI ,  Assistant  Chief 
Fuels  and  Lubrication  Division 
Aero  Propulsion  Laboratory 


"If  your  address  has  changed,  if  you  wish  to  be  removed  from  our 
mailing  list,  or  if  tne  addressee  is  no  longer  employed  oy  your 
organization  please  notify  AFWAL/POSF  ,  W-PAFB,  OH  Adijj  to  help 
us  maintain  a  current  mailing  list". 

Copies  of  this  report  should  not  be  returned  unless  return  is 
required  oy  security  considerations,  contractual  ou 1 1 o a t 1 ons ,  or 
notice  on  a  specific  document. 


ARTHUR  V.  CHURCHILL 
Chief,  Fuels  Branch 
Fuels  and  Lubrication  Division 
Aero  Propulsion  Laboratory 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

3EEORE  COMPLETING  .-CRM 

'  3L3QR  T  suMfiCS 

AFWAL-TS-M-2344  ,  Part  II 

2.  GOV’"  ACCSSiiON  NO. 

i^jec.R' j.m r's  h»jej 

a 

Laser  Diagnostic  Development  ar.d  Measurement 
and  Modeling  of  Turbulent  Tlowrields  of  Jets 
and  'v'akas ,  Part  IT 


7  AU 

1.  Krishnaxurthy , 

5  . 

0. 

Park,  ■; 

D.  J.  Wanrer,  and 

H. 

s . 

Cochran  j 

•  RSRRoamin a  3*ga*iirt;on  ume  *no  ascrsss 

University  of  Day  tor. 

Research  Institute 


3<,CGR  AM  £•_  -M 
a*Ea  i  «0*<  JNI*  NoMSE^S 

?.?.  i  52233?, 

Pro: .  =  3  3  43  , 


t  • .  Zr.H-ic. _ n  c  2  c  fflCc  n  AM  £  AMO  AOCR£SS  I  ’2.  : 

Aero  Propulsion  Laboratory  (AJWAL/'POS?)  ~~~ 

Air  Force  bright  Aeronautical  Laboratories  (A73C) 
•vrignt-Pattarson  Air  Force  Base.  Dhio  --- 


I  ’2.  P£?OB’  3aT£ 

|  June  19  3  3 

3.  nJm9£R  2*  3  *G  £5 


wA-aa.  ot  '-'-I  '*o<3*" 


.  :.N»a  .Report) 

>iic  Release,  Distribution  Unlimited 


'7  3lSTpi  9u  T?on  STaTSM£n"  'of  :/tm  «0«trac;  ,r»  3lock  20.  f  sitlmtfu  from  R  opart) 


'9.  <2V  *OROS  fConttnum  on  -ovoroo  j/cr#  r/  ‘tacfuarv  *rrj  :d*nttfy  ov  biocx  nuxno«r; 

Recirculating  Flowfields  Confined  Turbulent 

Isothermal  Flows  Cof lowing  Jets 

31uff-Bcdv  Near  Wake  Numerical  Predicti; 


20.  AflSTPAC"”  "Zontinum  on  **v«n«  ji<j#  ./  i«c«sierr  ertg  sv  3/ocx  itunflef) 

Finite-difference  numerical  computations  are  presented  to  pre- 
diot  the  isothermal,  turbulent,  recirculating  flowfields  in  a 
osnterbccy  combustor  configuration  which  involves  confined  dual 
coaxial  jet  mixing  in  the  near-wake  region  of  an  axi symmetric  | 

bluff  body.  The  calculations  based  upon  the  Reynclds-a-eraged  J 

Naviar-3 tokas  ecuatior.s  and  the  k-e  turbulence  model  consider  the  j 
influence  of  the  annular  and  central  flow  rates  on  the  nature  of 
the  flowfield  downstream,  of  the  bluff  bccv. 


DO  ’  1473  SCI-3N  3R  ■  MOV  IS  S  3SSOU4-S  JS„; 

5  c  C  L  3 1  - v  CL  AS3I  ff1  C  A  *'CN  2  c  i  S  3  Snt trod) 


LI 3 


SCCUBITV  CLASSIFICATION  OF  THIS  P AGCl'**T»«n  Omit  Enfnd) 


20.  continued 

The  isothermal  modelina  calculations  reported  herein  are  per¬ 
formed  with  the  *Teachina  Elliptic  Axisvmmetrical  Characteristics 
Heuristicallv*  (TEACH)  Code.  Present  calculations  have  employed 
the  standard  features  of  TEACH-type  numerics;  these  include  the  use 
of  primitive  variables  (velocity  components  and  pressure)  instead 
of  the  scream  f unction-vorticity  approach,  "hybrid"  upwind  dif¬ 
ferencing,  eddv-viscositv  approach  based  upon  the  k-e  model,  the 
Semi-Implicit  Method  for  Pressure-Linked  Equations  (SIMPLE) 
algorithm  for  the  pressure  field,  line-by-line  relaxation  and  tri- 
diaaonal  matrix  algorithm.  Additional  features  that  highlight  the 
present  calculations  are  the  power-law  differencing  scheme  which 
retains  the  diffusive  effects  for  a  larger  cell-Peclet  number 
range,  viz.,  -10  <_  Pe  <_  10,  than  the  "hybrid"  upwind  scheme;  and 
the  streamline-curvature  correction  in  the  k-e  model  which  involves 
a  curvature-dependent  (and  thus  nonconstant)  c-  (a  parameter  in  the 
"standard"  k-e  model  having  a  value  of  0.09) .  As  part  of  the 
numerical  investigations  of  the  centerbody  combustor  flowfields, 
the  predictions  of  the  modeling  without  and  with  these  additional 
features  are  compared. 

Present  -'^numerical  results  show  the  influence  of  the  annular  and 
central  jet  flow  rates  on  the  distributions  of  the  mean  and  rms 
velocity  fields  and  'the'- centerline  locations  of  stagnation  points. 
Moreover the  sensitivity  of  the  predicted  results  to  several 
aspects  of  the  modelina  is  considered.  These  include  the  dif¬ 
ferencing  schemes,  inlet  turbulence  lenath  scales,  streamline 
curvature  and  k-e  model  parameters. 

-  -.-The  predicted  results  demonstrate  the  complex  .nature  'of  the 
flowfield  interactions  in  the  near-wake  region  and  refine  the 
understanding  of  the  centerbody  combustor  flowfields.  The 
character  of  the  recirculating  flowfield  emergina  from  the 
numerical  predictions  when  the  near-wake  is  dominated  by  the 
annular  jet  is  in  conformity  with  the  experimental  observations.'-*' 
In  the  reverse-flow  reaion  behind  the  bluff  body,  the  present 
numerical  results  show  aood  agreement  with  the  recent  velocity  and 
concentration  measurements  and  with  the  annular- jet  data  in  the 
literature.  Also,  the  results  demonstrate  that  despite  the  com¬ 
plexity  of  the  centerbody  combustor  flowfield,  it  belonas  to  a 
wider  class  of  recirculatina  turbulent  flowfields  which  obey 
certain  similarity  considerations  for  the  mean  axial  velocity. 


Unclassified 


Cl»iTy  CLASSIFICATION  of  °  IGZ'Whan  En(ertd) 


1 


i 

N 

K 

[■ 

I  • 


- 

f 

y 

K. 

j 

I 

k 

* 

>/ 

I 


l 


* 

<• 


preface 


This  final  report  was  submitted  by  the  "'livers it/  of  Tw'''. 
jnler  Tor  tract  No.  ?3  36 13-73-C-200  5 .  The  project  was  s-ons  vre  d 
by  the  Air  -orce  Wright  Aeronautical  laoora tor ies ,  Aero 
Propulsion  Laboratory,  Wright-Patterson  Air  -orce  Base,  Ohio, 
under  Project  No.  3043,  Task  05,  Work  Unit  95.  Or.  William 
Roqueoore,  AFWAL/POSF,  was  Project  Engineer.  ""he  orogram  was 
managed  by  Dr.  Eugene  H.  Gerber  of  the  University.  Or.  L. 

Kr ishnanurthv  was  Principal  Investinator  of  the  analysis  task 
reported  in  this  volume.  This  research  task  was  initiated  in 
October  1930  and  completed  in  June  1932.  This  report  was  writte 
by  Dr.  L.  Kr ishnanurthv. 

The  Principal  Investigator  expresses  his  aooreciation  to: 

-  Dr.  W.  M.  Roquemore,  for  fruitful  discussions  during  this 
research  ; 

-  Dr.  R.  R.  Craig,  AFNAL/PORT,  for  making  available  a  ver¬ 
sion  of  the  TEACH-T  computer  program  and  for  helpful 

d iscuss  ions . 


Dr.  S.  0.  Park,  Mr.  D.  J.  Wahrer,  and  Mr.  H.  S.  Cochran, 
UDRI,  for  their  technical  contributions  to  the  numerical 
modeling  and  computing  activities; 


Ms.  A.  Mite,  Ms.  L. 
paring  this  report, 
n  i  c  a  1  ed  i  t  i  n  a  . 


Knox,  and  Mr.  M.  Jones,  UDRI,  for  ore 
and  Ms.  A.  Cochran,  UDRI,  for  tech- 


.  WL.  -  L  . _ ^ 


ABLE  OF  CONTENTS 


3 


3 


SECTION 


PAG 


i 


1  Introduction  i_ 

1.1  Background  i 

1.2  Obiectives  3 

1.3  Scope  of  Present  Work  5 

1.4  Outline  of  Report  7 

2  Governing  Equations  and  Solution  Procedure  9 

2.1  Differential  Equations  q 

2.2  Turbulence  Model  13 

2.3  Correction  for  Streamline  Curvature 
Effect 

2.4  Finite-Difference  Equations  19 

2.4.1  Exponential  Scheme  2S 

2.4.2  "Hybrid"  Upwind  Scheme  26 

2.4.3  Power-Law  Differencing  Scheme  23 

2.5  Boundary  Conditions  2^ 

.  2.6  Solution  Procedure  34 

2.6.1  Underrelaxation  35 

2.6.2  Convergence  36 

3  Preliminary  Numerical  Computations  37 

3.1  APL  Configuration  37 

3.1.1  Variations  of  Annular  and  Central 

Jet  Flow  Rates  39 

3.1.2  Sensitivity  Tests  42 

3. 1.2.1  Grid  Independence  of 

Solutions  42 

3. 1.2. 2  Influence  of  Inlet 

Velocity  Profile  44 

•  3. 1.2.3  Influence  of  Inlet 

Turbulence  Parameters  44 

3. 1.2. 4  Influence  of  the  Exit 

Boundary  Location  47 

3.2  UCI  Configuration  49 

3.2.1  Variations  of  Annular  and  Central 

Jet  Flow  Rates  49 

3.2.2  Sensitivity  Tests  50 

3. 2. 2.1  Influence  of  ae  50 

3.2.2. 2  Influence  of  cu  51 

3.2.2. 3  Other  Effects  53 

3.3  Similarity  Considerations  53 

3.3.1  Centerbody  Combustor 

Configurations  54 

3. 3. 1.1  UCI  Configuration  55 

3.3. 1.2  APL  Configuration  65 

3.3.2  Other  Flowfield  Configurations  65 

3.3.3  Species  Concentration  Fields  "’5 


v 


PREVIOUS  PAGE 
IS  BLANK 


TABLE  OF  CONTENTS  (Cont'i) 


SECTION  PAGE 

4  Refinements  and  Selected  Results  "9 

4.1  Influence  of  Differencing  Schemes  ~9 

4.2  COg  Mass  Conservation  30 

4.3  Influence  of  Inlet  Turbulence  Length 

Scales  Revisited  36 

4.4  Streamline  Curvature  Effects  in  Turbulence 

Modeling  93 

4.5  Comparison  of  Predictions  with  the  Newer 

Experimental  Results  98 

4.5.1  Centerline  Variation  of  the  Mean 

and  rms  Axial  Velocity  Fields  99 

4.5.2  Centerline  Variation  of  CO2 

Concentration  104 

4.5.3  Comparison  of  the  Measured  and 

Predicted  Streamlines  10” 

5  Conclusions  and  Recommer.dat  ions  113 

5.1  Conclusions  113 

5.2  Recommendations  for  Further  Activity  116 


References 


117 


LIST  OF  ILLUSTRATIONS 


IGURE 

PAGE 

1 

Computational  Oomain  and  Grid  Point 

Distribution  (Grid  A) 

20 

~> 

Schematic  of  Computational  Cell  Structure 

22 

3 

Influence  of  cu  on  the  Centerline  Staqnation 
Points  (UCI  Combustor) 

52 

4 

(a)  -  (g)  Mean  Axial  Velocity  Predictions 
for  the  (UCI)  Centerbody  Configuration 

56-59 

5 

(a)  -  (f)  Mean  Axial  Velocity  Predictions 
for  the  VUCI)  Centerbody  Conf iguration 

62-64 

6 

(a)  -  (c)  Mean  Axial  Velocity  Predictions' 
for  the  (APL)  Centerbody  Configuration 

66-67 

n 

Mean  Axial  Velocity  Measurements4^  for  the 
(APL)  Centerbody  Conf iguration 

68 

9 

Mean  Axial  Velocity  Measurements42  for  the 
(Abramovich)  Centerbody  Configuration 

69 

9 

Mean  Axial  Velocity  Measurements44  for  the 
Unconfined  Annular  Jet  Around  a  Centerbody 

71 

10 

(a)  -  (b)  Mean  Axial  Velocity4^  Measurements 
for  the  Unconfined  Annular  Jet  Around  a  Disk 

73-74 

11 

A  Composite  of  all  the  Predicted  and  Measured 
Results 

76 

12 

Radial  and  Axial  Distributions  of  COg  Mass 
Fraction 

81 

13 

Axial  Variation  of  the  Computed  Mass  Flow 

Rates  of  Air  and  COg 

83 

14 

An  Arbitrary  Modification  of  the  Grid  with 
Increased  Spatial  Resolution  in  the  Annular 
Region  (GRID  3) 

3  5 

15 

The  Relative  Importance  of  Diffusion  to  the 
Total  COg  Mass  Flow  Rate  at  Different  Axial 
Locations 

3^ 

LIST  OF  ILLUSTRATIONS  {  ZontM) 

I  SURE  PA.l 

In  Centerl  ine  Forward  and  Re a c  Stagnation  Joints  Vi 

17  Influence  of  Inlet  Turbulence  Length  Scale  and 

Curvature  Correction  on  the  Centerline 
Stagnation  Points  92 

13  Distribution  of  Curvature  Corrected  cu  15 

19  Effect  of  Streamline  Curvature  Correction  on 

the  Centerline  Mean  Axial  Velocity  Profiles  97 

20  Centerline  Mean  and  ms  Axial  Velocity  Profiles 

for  2  kg/s  Air  Flow  and  Eero  COo  Flow  IV) 

21  Centerline  Mean  and  ms  Axial  Velocitv  Profiles 

for  2  kg/s  Air  Flow  and  5  kg /hr  CD p  Flow  111 

22  Centerline  Mean  and  ms  Axial  Velocity  Profiles 

for  2  kg/s  Air  Flow  and  In  kg/hr  COg  Flow  102 

23  Centerline  COg  Mole  Fraction  for  2  kg/s  Air 

Flow  and  6  kg/hr  CO2  Flow  105 

24  Centerline  CO 9  Mole  Fraction  for  2  kn/s  Air 

Flow  and  1*5  kg/hr  CO2  Flow  196 

25  Streamline  Contours  for  2  kg/s  Air  Flow  and 

Zero  CO2  Flow  103 


LIST  OF  TA3I ES 


TABLE  PAG 

1  'I'j'Tr'! ing  Equations 

2  Boundary  Conditions  1 

3  Specification  of  e  at  the  Inlet  33 

4  Summary  of  Preliminary  Computational  Case 

Studies  3 0 

5  Influence  of  Grid  Spacings  on  (APLl  Centerline 

Stagnation  Distances  -13 

6  Influence  of  Inlet  Velocity  Profiles  on  (\?L) 

Centerline  Stagnation  Distances  45 

Influence  of  Inlet  Turbulent  Enemies  and  Length. 
Scales  on  (APL)  Centerline  Stagnation  Distances  4^ 

3  Influence  of  A  LAM  DA  and  TURBId  for  <3  kg/hr  CO n 

Flow  "  43 


LIST  OF  SYMBOLS 


a  Coefficient  for  the  conbi.oed  convective  and  diffusive 

flux  for  the  transport  of  any  variable 

b  Coefficients  arising  in  the  linearized  form  of  the 

source  term  see  Eq.  (14)  j 

c  Coefficients  arising  in  the  linearized  form  of  the 

source  tern  ,  see  Eq.  (14) 

cu  Constant  equal  to  0.09  in  the  "standard"  k- s  mod 

C  Convective  mass  flux  (=  pW  in  one-dinensional  or  .cm) 

Cj_  Constant  equal  to  1.44  in  the  k- e  model 

C a  Constant  equal  to  1.92  in  the  e  model 

D  Diffusive  mass  flux  (  =  ?/ 6z  in  one-dimensional  problem) 

centerbody  diameter 

f  Mixture  fraction 

H  Stagnation  enthalpy 

J  Combined  convective  and  diffusive  flux  see  Eq.  (19)  : 

k  Turbulent  kinetic  energy  Lsee  Eq.  (5)  : 

Kj_  Constants  appearing  in  the  expression  for  curvature- 
corrected  (nonconstant)  cy  .see  Eq.  (13)  j 

K2  Constants  appearing  in  the  expression  for  curvature- 
corrected  (nonconstant)  c u  ,see  Eq.  (13) 

Turbulence  length  scale  as  per  Eq.  (7b) 

l2  Turbulence  length  scale  as  per  Eq.  (9b) 

mair  Integrated  mass  flow  rate  of  air  per  unit  radian 

m^o 2  Inte<?rater5  mass  flow  rate  of  CC>2  per  unit  radian 

p  Static  pressure 


x 


r  „  Effective  exchange  coefficient  for  the  tame  lent  tr  ip.  or 
at  the  variable  >  _  see  Eq .  (  .1 '  ’ 

5  Characteristic  reference  lenc th  see  Tael-*  ■ 

3  p  Kr  "meeker  lelta  see  Eq  .  ;  10  > 

5z  Axial  -jrid  interval  between  ad  7  a  cent  nodal  points 

it  Normalized  radial  coordinate  3ue  to  Abr  anov  ich''  - 
(see  Paragraph  3.3.1) 

AW  Normalized  mean  axial  ve  1  oc  i  r  y  iue  to  An  r  a  mo  v  1  ;  h.  4  '  .  s  e  e 
Paraqar aph  3.3.1) 

s  Dissipation  rate  of  turbulent  kinetic  eneny  '  see  Eg. 

3  Azimuthal  coordinate  in  the  cylindrical  polar  aeometry 

\  Inlet  turbulence  lennth  scale  parameter  appear ino  in  E.i. 

(24)  (FORTRAN  variable  A  LAM.  DA  ) 

u  Absolute  viscosity,  molecular  viscosity  when  not 

subscripted 

v  Kinematic  viscosity 

0  Mass  density 

cjp  Effective  Prand 1 1/Schm id t  number  for  variable  >  .see  Table 
1  and  Eq  .  (  2  )  ] 

tw  Wall  shear  stress  (arises  in  wall-function  formulation  for 
the  specification  of  the  boundary  conditions  in  Paragraon 
2.4) 

>  General  dependent  variable  [see  Eq .  (1)[ 

SUBSCRIPTS 

A  Annular  jet  exit  plane,  air 

e  East  face  of  the  cell  control  volume 

E  East  nodal  point 

eff  Effective  value 


Eentral-jet  exit  plane 


f  ue  1 


i  Index  for  axial  computational  nodes,  components  in 

orthogonal  coordinate  directions 

in  Inlet  plane  of  the  computational  domain 

j  Index  for  radial  computational  nodes,  components  in 

orthagonal  coordinate  directions 

n  North  face  of  the  cell  control  volume,  coordinate  normal  to 

the  streamline 

N  North  nodal  point 

0  Oxidant 

P  General  computational  nodal  point 

s  South  face  of  the  cell  control  volume,  coordinate  along  the 

streamline 

S  South  nodal  point 

t  Turbulent 

w  West  face  of  the  cell  control  volume 

* 

W  West  nodal  point 

<s  General  dependent  variable 

SUPERSCRIPTS 

-  Time-Mean  value  (Reynolds  Averaging) 

'  Fluctuating  component  (Reynolds  Averaging) 

*  Dependent  variables  to  be  pressure-corrected  in  the  SIMPLE 

procedure  (see  Paragraph  2.6) 


SECTION  I 
INTRODUCTION 


This  three-oart  final  report  documents  the  research,  proa  ran 
performed  for  the  Air  Force  Wright  Aeronautical  laboratories, 

Aero  Propulsion  Laboratory,  by  the  University  of  Dayton.  The 
research  had  two  overall  objectives:  (a)  providing  profile  data 
that  can  be  used  to  evaluate  combustor  and  fuel  combustion  models 
and  (b)  evaluating  the  performance  of  combustor  models  and  dif¬ 
ferent  diagnostic  techniques  in  various  combustion  environments. 

The  technical  efforts  dealing  with  the  design  and  develop¬ 
ment  of  a  two-dimensional  laser  Doppler  anemometer  and  the 
exoerimental  data  collected  are  described  in  Part  I.  The  analy¬ 
sis  and  modeling  tasks  invoivinq  the  numerical  flowfield  predic¬ 
tions  and  their  comparisons  with  the  experimental  data  are 
described  in  this  volume,  Part  II.  Part  III  describes  the 
design,  development  and  performance  of  a  two-channel  time- 
resolved  laser  Raman  spectroscopy  system. 

1.1  BACKGROUND 

The  need  to  predict  complex,  rec irculating  turbulent 

flowfields  in  combustors  under  both  nonreacting  and  reacting  flow 

« 

conditions  has  provided  a  strong  incentive  for  the  develooment  of 
mathematical  models  and  numerical  procedures.  These  involve 
finite-difference  calculations  of  the  time-averaged  Navier-St okes 
equations  and  furnish  the  steady-state  predictions  of  the 
flowfields.  The  predictive  modeling  activities  have  made  use  of 
several  computer  programs  such  as  the  "Field  Relaxation  Elliptic 
Procedure"  (FREP)  code^  and  the  "Teaching  F.lliptic  Axisymmetrical 
Characteristics  Heur istically"  (TEACH)  code-.  Many  of  these 
modeling  investigations  have  be>  n  directed  at  the  centeroody 
combustor  in  operation  at  the  Air  Force  Wright  Aeronautical 
Laoor a  tor ies ,  Aero  Propulsion  Laboratory  ( AFWAL/PO 1  . 3“ 3  The 
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availability  of  centerbodv  combustor  experimental  data  from 
ongoing  research  programs,  involving  both  intrusive  and 
non  intrusive  diagnostic  techniques,  has  greatly  facilitated  the 
evaluation  and  validation  of  these  computer  codes. 

The  Aero  Propulsion  Laboratory  (APL)  ce.ntarbody  combustor 
configuration  represents  confined  dual  coaxial  jet  mixing  in  the 
near-wake  region  downstream  of  a  cylindrical  bluff  body. 

However,  the  internet  separation  is  much  larger  than  that 
encountered  in  tvoical  coaxial  jet  mixing,  since  the  con¬ 
figuration  involves  an  annular  jet  (flowing  between  the  confining 
outer  duct  of  0.254  m  diameter  and  the  centerbodv  of  0.14  n 
diameter)  and  a  central  jet  of  4.3-x-10"3  diameter.  This  wide 
a  separation  between  the  aets  and  the  concomitant  presence  of  the 
bluff-bodv  wake  in  the  mixing  region  have  not  been  encountered 
previously  in  numerical  modeling.  Thus,  the  APL  configuration 
provides  a  stringent  test  to  evaluate  the  predictive  capability 
of  several  candidate  modeling  codes. 

One  such  code  that  received  considerable  scrutiny  with 
respect  to  the  centerbody  flowfield  was  the  FREP  code,  which  was 
evaluated  in  earlier  modeling  research  studies  sponsored  by 
AFWAL/PO.3-5  one  of  these  previous  studies  clearly  established 
the  complex  nature  of  the  near-wake  flowfield  interactions  under 
different  annular  and  central  flow  rates.  5  The  character  of  the 
centerbody  flowfield  emerging  from  the  isothermal  predictions 
when  the  near  wake  was  dominated  by  the  annular  (air)  jet  or  the 
central  (COg)  jet  was  found  to  conform  to  the  APL  experimental 
observations  and  with  the  heuristic  flowfield  descriotions 
suggested  therefrom.'5  Although  the  predicted  results 
demonstrated  the  capability  of  the  FREP  code  to  provide 
qualitatively  correct  trends  in  the  flowfield  behavior,  a 
comparison  of  the  predictions  with  the  measured  velocity^  and 
the  concentration^  data  revealed  only  fair  to  poor  quantitative 
agreement.  The  major  shortcoming  of  the  FREP  modeling  which 
contributed  to  the  poor  agreement  between  the  predicted  and 
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measured  results  was  known  to  be  in  the  turbulence  model.  The 
numerical  calculations 5  had  employed  the  obviously  crude  model 
which  assumed  a  constant  value  of  turbulent  eddy  viscosity,  since 
the  FREP  code  had  encountered  serious  numerical  difficulties  with 
the  well  known  two-equation  (k-e )  turbulence  model. 

The  challenge  inherent  in  the  centerbody  configuration  and 
the  availability  of  APL  diagnostic  data^-2  thereon  evoked  the 
interest  of  a  number  of  modelers  (e.g.,  at  AIRESEARCH , 13 
MASA-Lewis^  and  Pratt  and  Whitney^'^)  in  examining  the  center- 
body  configuration  with  their  computer  codes.  These  codes  are 
based  upon  or  variants  of  the  TEACH  program,  originally  developed 
by  Gosman  and  Pun 2  and  subsequently  tested  and  validated  by  a 
number  of  investigators .  15-10  Therefore,  the  TEACH  code  appeared 
to  be  an  attractive  candidate  for  use  in  our  model  validation 
studies  with  respect  to  the  centerbody  combustor  flowfields.  A 
version  of  the  computer  code,  TEACH-T,  was  made  available  to  us 
by  the  Aero  Propulsion  Laboratory,  Ramjet  Technology  Branch 
(AFWAL/PORT) . 20  Since  this  code  contained  an  operational  k-e 
turbulence  model,  it  became  of  interest  to  test  this  code  for  the 
centerbody  configuration. 

1.2  OBJECTIVES 

The  analysis  task  was  concerned  with  providing  the  necessary 
theoretical  framework  for  the  interpretation  of  the  diagnostic 
data  on  the  velocity  and  concentration  fields  made  available  from 
the  ongoing  APL  experimental  program.  This  required  numerical 
predictions  of  the  flowfield  variables  (for  both  mean  and 
fluctuating  components)  for  various  experimental  conditions  and 
their  comparison  with  the  experimental  data.  This  procedure  was 
also  expected  to  facilitate  the  evaluation  of  the  performance  of 
combustor  models  and  diagnostic  techniques.  The  availability  of 
computer  programs  that  are  reasonably  successful  in  simulating 
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the  desired  flowfields  was  essential  to  the  completion  o£  the 
analysis  task. 

It  was  recognized  at  the  outset  that  the  EP.E?  code  was 
deficient  in  terms  of  the  turbulence  model  and  that  recourse  muse 
be  made  to  other  finite-difference  codes  availaole  in  the  public 
domain.  The  availability  of  the  TEACH-T  program  from  AFWAL/PORT 
essentially  determined  the  course  of  the  research  program. 
Accordingly,  the  implementation  of  the  analysis  task  was 
considered  in  terms  of  the  following  specific  -objectives: 

a.  Modify  the  AFWAL 7 PORT  version  of  the  TEACH-T  program  for 
the  centerbodv  combustor  flowfield,  make  it  operational,  and 
verify  its  capability  to  predict  the  isothermal  flowfields. 

b.  Carry  out  performance  assessments  of  the  code  through 
parameter  optimizations,  modif ications  and  improvements,  in 
comparison  with  other  model  predictions  as  well  as  experimental 
data . 

c.  Provide  the  numerical  predictions  of  the  velocity  and 
concentration  fields  for  the  desired  experimental  conditions. 

The  successful  accomplishment  of  these  objectives 
necessitated  a  close  and  continuous  involvement  with  the  APL 
experimental  program.  It  is  easy  to  recognize  the  need  for  such 
interaction  for  the  comparison  of  model  predictions  with 
experimental  data.  In  addition,  the  numerical  modeling  required 
input  parameters  appropriate  for  the  desired  experimental 
conditions.  The  selection  of  these  input  requirements  (in  terms 
of  the  geometric,  fluid  dynamic,  and  chemical  parameters,  and  of 
the  boundary  conditions  involving  velocity,  temperature, 
pressure,  and  species  mass  fractions)  would  have  to  be  based  upon 
data  made  available  from  the  experimental  program.  Thus, 
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fulfillment  of  the  above  objectives  depended  on  the  timely 
availability  of  appropriate  experimental  data. 

1.3  SCOPE  OF  PRESENT  WORK 

After  a  successful  modification  of  the  AFWAL/PORT  version  of 
the  TEACH-T  code  for  the  centerbody  combustor  flowfield  early  in 
the  program,  preliminary  computations  of  the  isothermal  mixing  of 
2  'Kg/s  annular  air  flow  and  2.22xl0“3  kg/s  (8  kg/hr)  central  CO2 
flow  revealed  numerous  numerical  instabilities,  convergence 
failures,  and  nonphysical  flowfield  features.  Since  the  code  was 
written  explicitly  for  the  reacting-flow  computations  and  also 
required  the  correct  specifications  of  the  excess  air  ratio  to 
determine  the  air  flow  rate  from  the  specified  fuel  flow  rate, 
substantial  changes  were  needed  to  employ  the  code  for 
nonreacting  calculations.  Therefore,  the  initial  phase  of  the 
modeling  research  was  devoted  to  making  the  code  operational  for 
the  isothermal  predictions  of  the  centerbody  flowfields. 

When  the  required  changes  were  implemented  in  the  code,  the 
numerical  computations  were  found  to  proceed  smoothly,  resulting 
in  physically  acceptable  solutions.  Converged  results  of  the 
calculations  for  2  kg/s  air  flow  and  several  CO2  flows  (such  as 
1.11x10-3,  2.22xl0“3f  3.33x10-3,  4.44x10-3,  and  5xl0~3  kg/s)  were 
obtained  and  compared  with  the  available  experimental  data  on  the 
velocity  field. 30  In  addition,  calculations  were  made  for  the 
extremely  small  CO2  flow  rate  of  2.8xlQ-7  kg/s  (lxl0“3  kg/hr)  and 
the  results  were  considered  to  be  a  good  approximation  of  those 
of  the  zero  CO2  flow  rate  case.  (This  approach  was  necessary 
since  the  initial  programming  of  the  code  did  not  allow  the 
specification  of  zero  CO2  flow  rate.)  These  results  were 
compared  with  the  experimental  data  for  zero  CO2  flow  rate. 

These  comparisons  showed  good  qualitative  agreement  for  the  mean 
axial  velocity  and  turbulence  intensity  profiles  along  the 
centerline.  The  quantitative  agreement  for  the  centerline 
locations  of  the  stagnation  points,  however,  was  inadequate. 


Investigation  of  the  causes  of  this  discrepancy  required  a  series 
of  sensitivity  tests  and  optimization,  since  the  modeling 
depended  on  a  number  of  factors.  These  were  the  spatial 
resolution  of  the  computational  grid,  the  specifications  of  the 
inlet  velocity  profile,  turbulence  intensity  and  length  scale, 
the  prescription  of  the  values  of  the  constants  in  the  turbulence 
model,  and  the  specification  of  the  location  of  the  exit  boundary 
of  the  computational  domain.  The  last  item  is  always  arbitrary 
and  subject  to  a  trade  off  between  the  available  computer 
resource  and  desirable  solution  accuracy.  Therefore,  the  next 
phase  of  the  modeling  research  examined  the  questions  of 
parameter  sensitivity  and  optimization. 

Concurrent  with  the  experimental  program  on  the  APL  center- 
body  combustor,  investigations  have  been  underway  on  a  centerbody 
combustor  of  roughly  1/5-scale  model  of  the  APL  configuration  at 
the  Combustion  Laboratory  of  the  University  of  California,  Irvine 
(UCI). 21,22  since  experimental  data  were  becoming  available  for 
this  combustor,  it  was  of  interest  to  obtain  the  numerical 
predictions  thereon  with  the  TEACH  code.  It  was  expected  that 
such  calculations  would  provide  information  on  the  question  of 
combustor  scaling  (as  between  the  large-scale  APL  and  the  small- 
scale  UCI  configurations).  Moreover,  these  calculations  also 
served  to  examine  additional  sensitivity  tests  and  the  question 
of  similarity7  of  the  radial  distribution  of  the  mean  axial 
velocity.  Accordingly,  these  aspects  represented  another  phase 
of  the  modeling  research. 

The  sensitivity  influence  and  parameter  optimization  did  not 
appear  to  reduce  the  disparity  between  the  predicted  and 
measured^  centerline  stagnation  point  locations.  Moreover, 
while  the  model  calculations^  at  Pratt  and  Whitney  also  resulted 
in  the  underprediction  of  the  forward  stagnation  point  and  the 
overprediction  of  the  rear  stagnation  point,  our  results 
indicated  greater  degrees  of  under-  and  over-prediction.  The 
search  for  an  explanation  of  this  somewhat  disturbing  anomaly  led 


6 


to  a  reexamination  of  the  influence  of  the  inlet  length  scale  and 
also  to  the  question  of  incorporating  a  correction  in  the 
turbulence  model  to  account  for  streamline  curvature.  These 
aspects  formed  the  final  phase  of  the  modeling  research  within 
the  scope  of  the  present  program. 

The  several  phases  of  modeling  research  outlined  above 
conformed  to  the  first  two  specific  objectives  indicated  in 
Paragraph  1.2.  The  AFWAL/PORT  version  of  the  TEACH  code  may  now 
be  said  to  be  operating  satisfactorily,  with  suitable  parameter 
optimizations  and  improvements  having  been  implemented.  As  seen 
in  Section  4,  the  present  status  of  our  predictive  modeling  is 
good.  However,  the  modeling  research  remains  an  ongoing 
activity,  particularly  with  respect  to  the  third  objective. 
Extensive  experimental  data  have  become  available  recently  from 
the  UCI  combustor. 23  These  data  as  well  as  the  isothermal  LDA 
data  of  the  APL  combustor  (seen  in  Part  I  of  this  report) 
necessitate  further  computational  efforts  for  the  comparisons 
with  the  numerical  predictions.  These  additional  computations 
would  serve  to  examine  the  question  of  combustor  scaling  and  help 
elucidate  the  distinctions  between  the  small-  and  large-scale 
configurations . 

1.4  OUTLINE  OF  REPORT 

The  governing  equations  and  the  solution  procedure  of  the 
TEACH-T  program  are  discussed  briefly  in  Section  2.  Aspects  of 
the  two-equation  turbulence  model,  the  numerical  algorithms, 
convergence  criteria  etc.  are  highlighted  in  that  section.  In 
Section  3  the  preliminary  computational  case  studies  are 
identified  and  the  modeling  details  are  discussed.  The  results 
of  the  refined  computations  are  presented  and  discussed  in 
Section  4.  The  numerical  predictions  are  compared  first  with 
the  earlier  LDA  data.  10  These  comparisons  facilitated  the 
implementation  of  a  number  of  computational  refinements.  With 
the  availability  of  the  newer  results  documented  in  Part  I, 
additional  comparisons  thereof  with  the  computed  results  are 
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SECTION  2 

GOVERNING  EQUATIONS  AND  SOLUTION  PROCEDURE 


The  numerical  solution  of  a  system  of  partial  differential 
equations  describing  the  conservation  laws  of  mass,  momentum, 
energy,  and  chemical  species  forms  the  basis  of  the  predictive 
modeling  of  the  turbulent,  recirculating  flowfields  in  the  cen- 
terbody  combustor  configuration.  This  is  accomplished  in  the 
TEACH  code  by  a  finite-difference  computational  procedure  to 
solve  the  time-averaged  Navier-Stokes  equations.  A  brief 
discussion  of  the  theoretical  and  computational  aspects  of  the 
TEACH  code  is  given  below.  Further  details  are  available  in 
References  2  and  15. 

2.1  DIFFERENTIAL  EQUATIONS 

The  TEACH-T  computer  program  is  written  for  steady-state 
flowfields  in  two-dimensional  {planar  or  axisymmetric) 
geometries.  The  formulation  entails  an  elliptic  system  of 
equations  for  properly  describing  the  recirculating  flows.  The 
code  has  been  extended  by  a  number  of  investigators  to  apply  to 
three-dimensional  geometries  and  can  also  be  easily  modified  to 
include  parabolic  flows.  However,  only  the  axisymmetric  geometry 
and  the  elliptic  formulation  are  of  interest  in  the  present 
program.  Furthermore,  since  the  present  computations  are  con¬ 
cerned  only  with  isothermal  flowfields,  the  governing  equations 
do  not  take  the  chemical  source  terms  into  account. 

The  numerical  treatment  of  the  Navier-Stokes  equations  in 
the  TEACH  code  involves  the  primitive  (pressure  and  velocity) 
variables  instead  of  the  stream  function-vorticity  formulation 
employed  in  the  FREP  code.l  Therefore,  a  direct  solution  of  the 
velocity  and  pressure  fields  cannot  be  avoided.  The  code  uses  a 
special  procedure  called  the  Semi-Implicit  Method  for 
Pressure-Linked  Equations  (SIMPLE) 24  algorithm  to  solve 
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explicitly  for  the  velocity  and  pressure  fields.  Further  details 
concerning  the  underlying  theory  and  the  computational  procedure 
are  available  in  References  15/  24  and  25. 


For  the  ax isymmetr ic  configuration  of  interest  here,  we 
write  the  governing  differential  equations  in  the  cylindrical 
polar  coordinates.  In  view  of  the  notation  adopted  in  the  APL 
experimental  program,  our  present  nomenclature  for  the  coor¬ 
dinates  and  velocity  components  is  different  from  that  in  the 
original  TEACH  code  formulation . Thus  U,  V  and  W  are  the  time- 
mean  velocity  components  in  the  radial  (r),  azimuthal  (8),  and 
axial  (z)  directions.  For  the  axisymmetric  (in  the  mean) 
flowfield,  the  flow  variables  have  no  explicit  dependence  on  9 
and  for  the  case  of  no  swirl  considered  here,  V  is  zero.  The 
governing  equations  for  all  the  dependent  variables  can  be 
expressed  in  the  general  form 


1  fa_ 

r  3  r 


(prU*  )  +  (prW*  )  -  ^ 


(rr  11)  -  L_ 
'  <fr  3  r '  3  z 


3  4>  v  J  _ 
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where  <j>  denotes  any  dependent  variable  (time-mean  value).  In 
Eq.(l)  is  the  source  term  for  the  variable  ij>  which  includes 
true  source  terms  (such  as  those  due  to  chemical  reactions)  as 
well  as  the  terms  not  covered  by  the  first  four  terms 
(representing  the  convective  and  diffusive  contributions), 
is  the  effective  exchange  coefficient  for  the  transport  of  the 
variable  $  and  is  given  by 

r $  =  u ef f/o*  ( 2 ) 

where  ueff  is  the  effective  viscosity  in  the  flowfield  and  is 
the  appropriate  effective  Prandtl/Schmidt  number  for  each  . 
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In  a  general  formulation  of  the  problem  dealing  with  three- 
dimensional,  reacting  flowfields,  the  dependent  variable  }>  de¬ 
notes  U,  V,  W,  the  stagnation  enthalpy  H,  fuel  mass  fraction  Yp, 
oxidant  mass  fraction  Y q,  mixture  fraction  f,  turbulent  kinetic 
energy  k,  its  dissipation  rate  e,  the  mean  square  fluctuations  of 
concentrations  (such  as  f'2,  y.J.2  and  Y^2),  and  any  other  variable 
such  as  YpY^  (e.g.,  see  Reference  26).  Present  calculations, 
however,  deal  with  the  nonreacting  flowfields  (as  arising  in  the 
turbulent  mixing  of  the  annular  air  stream  and  the  central  CO2 
stream  in  the  centerbody  configuration).  Therefore,  for  the  spe¬ 
cies  conservation  equation  dealinq  with  the  variable  Yp  (which  is 
now  a  passive  scalar),  the  source  term  due  to  chemical  reaction 
vanishes.  For  the  energy  conservation  equation  in  the  AFWAL/PORT 
version  of  the  code,  the  dependent  variable  is  the  static 
enthalpy.  For  the  isothermal  flowfields  under  consideration 
here,  it  is  not  necessary  to  compute  the  static  enthalpy. or  the 
temperature.  Thus,  the  dependent  variables  of  present  interest 
are  U,  W,  p,  Yp,  k  and  e.  The  last  two  variables  are  discussed 
in  detail  in  Paragraph  2.2.  The  appropriate  values  of  and  the 
details  of  for  the  dependent  variables  of  interest  are  shown 
in  Table  1.  The  one  remaining  unknown  variable  is  the  pressure 
p,  for  which  there  is  no  additional  governing  equation.  The 
SIMPLE  algorithm2**  provides  a  special  procedure  for  obtaining  the 
pressure  field  in  the  TEACH  code.  Briefly,  the  momentum 
equations  are  first  solved  by  estimating  a  pressure  field,  and 
then  the  computed  velocity  field  is  used  to  correct  the  pressure 
field  by  ensuring  that  the  continuity  equation  is  satisfied. 

The  effective  viscosity  yeff,  appearing  in  Eq.  (2)  remains 
to  be  determined.  In  our  formulation,  ueff  is  the  exchange  coef¬ 
ficient  for  the  transport  of  momentum  in  the  turbulent  flowfield. 
Thus,  for  the  dependent  variables  U  and  W,  r ^  is  given  by  ueff, 
which  is  different  from  and  much  larger  than  the  molecular 
viscosity  u  of  the  fluid  of  interest.  Since  we  deal  with  the 
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TABLE  1 

GOVERNING  EQUATIONS 


time-averaged  equations,  the  process  of  time  averaging  introduces 
correlations  of  fluctuating  velocity  components  as  unknowns  in 
the  mean  momentum  equations.  Also,  the  numerical  solutions 
involve  computational  cell  sizes  that  are  much  larger  than  the 
length  scales  of  importance  in  the  actual  turbulent  motion.  The 
method  of  approximating  the  unknown  correlations  in  terms  of  the 
known  quantities  gives  rise  to  turbulence  modeling.  This  aspect 
is  discussed  next. 

2.2  TURBULENCE  MODEL 

The  effective  viscosity  ueff  appearing  in  Eq .  (2)  is 

expressed  as  the  sum  of  the  molecular  (or  the  laminar)  viscosity 
u  of  the  fluid  and  a  turbulent  eddy  viscosity  ut,  i.e., 

ueff=u+uf  (3) 

The  TEACH  code  employs  the  so-called  two-equation  turbulence 
model  which  relates  ut  to  two  scalar  properties  of  turbulence,  k 
and  e  as 

u  t  =  cu  o  k2/e  ,  ( 4 ) 

where  c^  is  usually  taken  to  be  a  constant  equal  to  0.09.  The 
above  procedure  involves  the  introduction  of  two  partial  dif¬ 
ferential  equations  for  k  and  e  which  are  solved  together  with 
the  conservation  equations  for  mass,  momentum,  and  energy.  We 
note  that  the  inclusion  of  these  additional  equations  has  been 
anticipated  in  the  general  formulation  of  Eq .  (1)  and  in  Table  1. 
It  is  stressed  here  that  the  pressure  fluctuations  are  not  con¬ 
sidered  and  that  the  source  terms  S^,  for  W,  u,  k,  and  e  in  Table 
1  are  those  appropriate  for  incompressible  flow,  since  the 
AFWAL/PORT  version  of  the  code  is  explicitly  written  for  that 
case.  Therefore,  overbars  for  p  and  p  are  superfluous.  The 
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required  modifications  to  compressible  flow  may  be  found  in 
References  26  and  27. 

Formally,  the  turbulence  kinetic  energy  k  is  taken  to 
characterize  the  intensity  of  turbulent  fluctuations ,  i.e., 

k  i  TTui/2  (5) 

where  Uj_'s  are  the  fluctuating  velocity  components  in  the  three 
orthogonal  directions.  In  Eq.  (5)  the  familiar  convention  of 
summation  on  repeated  indices  is  adopted.  A  number  of  authors 
have  modeled  the  transport  equation  for  k  (see  e.g.,  Launder  and 
Spalding^S).  The  modeled  equation  appropriate  for  the  axisym- 
metric  geometry  without  swirl  considered  here  is  shown  in  Sq. 

(  1 )  and  Table  1 . 

The  rate  of  viscous  dissipation  of  turbulence  kinetic  energy 
is  the  second  turbulence  scalar  of  interest.  For  large  Reynolds 
numbers,  this  variable  may  be  defined  as 
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3u .  3u  • 
3  3 


(6) 


where  v  is  the  kinematic  viscosity  (s  u/p),  the  Xj's  are  the 
distances  along  the  coordinates  r,  3  and  z,  and  the  summation 
convention  is  used.  The  modeled  form  of  the  transport  equation 
for  s  given  by  Eq.  (1)  is  that  appropriate  for  high-Reynolds- 
number,  ax i symmetric  recirculating  flows  without  swirl. 

The  standard  k-£  model  discussed  in  the  foregoing  enables 
the  calculation  of  the  turbulent  eddy  viscosity  through  Eq.  (4). 
However,  the  knowledge  of  k  and  £  alone  does  not  describe  the 
structure  of  the  turbulence.  Still  it  is  possible  to  deduce  some 
information  on  the  characteristic  length  and  time  scales  of 


turbulence  from  k  and  e.  Thi  procedure  is  facilitated  by  the 
use  of  the  Prandtl-Kolmogorov  relation, 


from  which  the  length  scale  Zj_  is  related  to  through  Eq.  (4)  as 

Zx  =  *3/2/ e.  (7b) 

The  corresponding  time  scale  t]_  which  is  given  by  Z^/k^/2  is 
obtained  from 


t!  =  k/  e.  (8) 

The  formulation  illustrated  by  Sqs .  (7a)  and  (7b)  has  found 
application  in  some  earlier  studies  (e.g.,  see  References  29  and 
30).  Often,  instead  of  Eqs.  (7),  an  alternative  formulation  is 
employed  (e.g.,  see  Reference  31),  according  to  which 


1/2 

ufc  -  p  K  (9a) 

and 

Z  =  c  k3/2/e  .  ( 9b ' 

2  u 

We  note  that  the  AFWAL/PORT  version  of  the  TEACH  Code 
employs  the  former  formulation,  while  the  recent  calculations  of 
Sturgess  and  Syed^  are  based  upon  the  latter  formulation.  An 
inspection  of  Eqs.  (7)  and  (9)  reveals  that  the  length  scales 
and  I2  differ  by  a  factor  of  cy.  This  distinction  must  be  taken 
into  account  when  the  inlet  boundary  condition  for  e  (see 
paragraph  2.4)  is  specified  through  the  prescription  of  the  inlet 
length  scale.  Indeed,  this  is  the  only  circumstance  when  the 
length  scale  comes  into  play  in  the  calculations,  since  ut  is 
calculated  by  Eq.  (4)  for  the  distributions  of  k  and  e  obtained 
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by  solving  Eq.  (1)  and  the  set  of  Eqs .  (7)  or  (9)  is  not  directly 

involved  in  the  computations. 

2.3  CORRECTION  FOR  STREAMLINE  CURVATURE  EFFECT 

It  is  well  known  that  the  streamline  curvature  has  a  strong 
influence  on  the  shear-flow  turbulence.  This  topic  has  been 
extensively  discussed  by  Bradshaw, 32  whose  review  of  some  experi¬ 
mental  studies  shows  that  the  turbulent  shear  stress  and  the 
degree  of  anisotropy  between  the  normal  stresses  are  very  sen¬ 
sitive  to  curvature.  Numerical  predictions  of  the  size  of  the 
recirculation  regions  encountered  in  such  turbulent  recirculating 
flows  as  the  centerbody  and  sudden  expansion  geometries  also 
appeared  to  depend  strongly  on  the  turbulence  acti/ity  in  the 
curved  shear  layers  bordering  the  recirculation  region.  The 
standard  k- s  model  discussed  in  Paragraph  2.2  does  not  account 
for  streamline  curvature  effects.  It  can  therefore  be  expected 
that  this  failure  may  be  partly  responsible  for  the  discrepancy 
between  the  predicted  and  measured  recirculation  lengths. 

Curvature  modifications  to  the  standard  k- e  model  have  been 
attempted  by  a  number  of  authors  (e.g.,  see  References  33  through 
35).  All  these  corrections  are  ad  hoc,  and  their  physical  basis 
and  range  of  applicability  are  open  to  doubt.  Further  research 
is  needed  to  provide  a  more  rigorous  framework  for  the  incor¬ 
poration  of  curvature  correction  in  the  turbulence  models.  In 
our  present  program,  however,  we  have  introduced  a  curvature- 
dependent  (and  hence  nonconstant)  c u  into  the  standard  k- e  model, 
along  the  .lines  suggested  by  Leschziner  and  Rodi.33  This  modifi¬ 
cation  is  briefly  described  below. 

The  streamline  curvature  modification  in  Reference  33  is 
based  upon  the  algebraic  stress  model  of  Gibson36  which  essen¬ 
tially  consists  of  a  set  of  algebraic  equations  resulting  from 
the  deletion  of  all  transport  terms  in  the  differential  equations 


16 


governing  the  transport  of  the  Reynolds  stresses. 37  These 
equations  can  be  shown  in  a  general  form  as 
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where  Pj_j  is  the  rate  of  production  of  the  turbulent  stress  uj_Uj  , 
(shown  in  Table  1)  is  the  generation  of  turbulence  energy  by 
the  interaction  of  mean  velocity  gradients  and  turbulence 
stresses,  6j_j  is  the  Kronecker  delta,  and  a  and  3  are  constants 
having  the  values  of  1.5  and  0.6  respectively  (as  in  Reference 
37).  With  the  assumption  of  local  equilibrium  of  turbulence 
energy  between  production  and  dissipation  (expressed  by  =  e) , 
Eq.  (10)  takes  the  simplified  form 


uiuj  =  1-3 
~ 1c  ae 
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(11) 


Equation  (11)  reveals  the  effect  of  curvature  on  the 
stresses  if  (i,j)  are  taken  as  the  streamline  coordinates  (s,n) 
where  s  is  the  coordinate  along  the  streamline  and  n  is  the 
coordinate  along  the  local  normal  to  the  streamline.  Then  Pj_j 
becomes 
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In  Eq.  (12)  Rc  is  the  local  radius  of  curvature,  r  is  the  radial 
distance  from  the  axis  of  symmetry  (for  the  axially  symmetric 
flow  considered  here)  and  Us ,  Un  and  Ur  are  respectively  the 
velocity  components  along  the  s,  n,  and  r  directions.  r,»'e  can 
determine  the  stresses  u-,  "uj  and  unus  from  Eqs.  (11)  and  (12)  if 
the  mean-velocity  field,  k,  and  e  are  known  or  determinable. 
However,  the  shear  stress  unus  is  the  dominant  stress  term  in  the 
momentum  equations  in  a  curved  shear  layer.  Therefore,  it  would 
suffice  to  focus  our  attention  on  usun ,  obtain  an  expression 
which  relates  this  term  to  the  respective  rate  of  strain  (  3(JS/  3n 
+  Us/Rc),  and  thereby  extract  a  curvature-dependent  c  u.  This 
result  is  shown  as 

2  v2  3U„  U  U 

Cu  5  -  K1K2^1+8K1  ^<T?T  +  TT>  TT  ,  <13> 

£  C  C 

where  Kj_  s  (1-3) /a  and  K2  =  (2/3)  (l-a-3)/a. 

Equation  (13)  shows  that  cu  is  no  longer  a  constant  but  is  a 
function  of  streamline  curvature.  If  the  terms  3US/ 3n  and  Ug/Rc 
are  known  along  with  k  and  s,  cu  is  evaluated  from  Eq.  (13).  In 
Reference  33  3US/ 3n  and  Us/Rc  are  evaluated  from  W,  U  and  their 
derivatives  in  z  and  r  directions,  as  part  of  the  main  solution 
algorithm.  These  details  are  not  included  here. 

Before  we  conclude  this  topic,  it  is  worthwhile  to  emphasize 
the  ad  hoc  nature  of  the  cu  correction  derived  here.  According 
to  the  specified  values  of  a  and  3,  K]_K2  =  -0.13.  Thus,  when  Eq . 
(13)  is  used  with  the  k- £  model  and  for  negligible  correction  for 
curvature  effects  Lthe  denominator  in  Eq.  (13)  approaches  unity  j, 
cu  reduces  to  0.13,  in  contrast  with  the  standard  value  of  0.09. 
To  avoid  this  inconsistency.  Reference  33  assumes  -K]_K2/c  u  =  1 
(instead  of  1.5)  in  the  actual  calculations.  But  the  value  of 
3K^  in  the  denominator  of  Eq.  (13)  is  determined  from  the  spe¬ 
cified  values  of  a  and  3.  This  inconsistency  clearly  introduces 
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a  more  con- 


certain  arbitrariness;  as  indicated  in  Reference  35, 
sistent  approach  would  be  to  determine  8K^  from  the  condition 
-K1K2  =  0.09.  Another  feature  in  Eq.  (13)  that  deserves  comment 
is  that  c.j  may  become  negative  when  (  3US/  3n )  /  ( Us/Rc )  is  negative 
and  sufficiently  large  in  magnitude.  This  will  be  physically 
absurd  in  view  of  Eq.  (4).  Thus,  the  actual  computational  algo¬ 
rithm  imposes  an  arbitrary  positive  lower  bound  on  c y.  This 
aspect  is  discussed  later  in  Paragraph  4.4. 


2.4  FINITE-DIFFERENCE  EQUATIONS 

The  differential  equations  describing  the  conservation  laws 
are  represented  by  Eq.  (1).  The  solution  of  these  equations  is 
obtained  numerically.  This  procedure  consists  of  specifying  a 
computational  grid  distribution  for  the  flow  domain,  obtaining  a 
set  of  finite-difference  algebraic  equations  derived  from  discre¬ 
tizing  the  differential  equations  on  all  the  grid  nodes  and 
solving  the  algebraic  equations  by  standard  numerical  methods. 

9 

The  finite-difference  equations  can  also  be  directly  obtained 
from  a  control-volume  analysis  of  the  conservation  laws  (see 
Reference  15  for  details).  The  accuracy  of  the  numerical 
solution  depends  on  how  closely  the  set  of  finite-difference 
equations  approximates  the  original  differential  equations.  In 
general  this  is  governed  by  the  number  of  computational  grid 
nodes  which  represent  the  flowfield. 

Figure  1  shows  the  computational  domain  and  the  grid-point 
distribution  initially  adopted  in  the  numerical  calculations. 

The  chosen  grid  consists  of  41  axial  nodes  and  34  radial  nodes, 
and  viewed  in  the  r-z  plane,  is  regular  and  rectangular  with 
nonuniform  grid  spacing  in  both  axial  and  radial  directions.  As 
pointed  out  earlier,  the  location  of  the  exit  boundary  is 
aroitrary  and  the  present  location  of  30  cm  from  the  face  of  the 
centerbody  has  been  based  upon  extensive  computational 
experience.  To  ensure  adequate  spatial  resolution  in  flowfield 
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Figure  1.  Computational  Domain  and  Grid  Point  Distribute 
(GRID  A) . 


regions  with  large  gradients  of  the  flow  variables,  the  grid- 
point  distribution  is  more  densely  populated  toward  the  center- 
body  face  and  the  axis  of  symmetry  than  toward  the  exit  boundary 
and  the  duct  wall.  All  the  salient  locations  such  as  the  axis  of 
symmetry,  the  central-tube  boundary,  the  centerbody  tou  boundary 
and  trailing  face,  duct  wall  boundary,  and  the  exit-plane 
boundary  are  located  midway  between  the  adjacent  grid  nodes. 

Figure  2  shows  the  schematic  of  the  grid  arrangement 
employed  in  the  formulation  of  the  finite-difference  equations. 

All  the  dependent  variables,  except  V  and  W,  are  referred  to  at 
the  grid  nodes  (e.g.,  the  node  denoted  by  P) .  0  and  W  are  calcu¬ 

lated  at  locations  midway  between  the  grid  nodes,  as  illustrated 
by  the  dotted  lines.  In  other  words,  a  staggered -grid  arrange¬ 
ment  has  been  employed.  This  has  the  advantages  of  (a)  eval¬ 
uating  the  pressure  gradients  (which  drive  the  velocities  0  and 
W)  at  the  locations  where  the  velocities  are  calculated  and  (b) 
determining  the  velocities  where  they  are  needed  for  the  calcula¬ 
tion  of  the  convective  fluxes  into  and  out  of  each  cell  (viz.,  at ' 
the  cell  boundaries)  . 

The  set  of  finite-difference  equations  for  Eq.  (1), 
excluding  the  continuity  equation  which  receives  special  treat¬ 
ment,  is  derived  by  integrating  Eq.  (1)  over  a  control  volume 
surrounding  each  grid  node.  This  procedure  requires  appropriate 
assumptions  to  describe  the  relation  between  the  nodal  value  ?p , 
the  rate  of  creation/destruction  of  $  within  the  control  volume, 
and  the  transport  of  convective  and  diffusive  fluxes  of  $  across 
the  boundaries  of  the  control  volume.  The  source  term  S ^  in  Eq. 
(1)  (see  Table  1  for  the  form  of  S  ^  for  various  dependent 
variables)  can  be  represented  in  a  linearized  form  as 

/  S  dV-  =  b  a  +  c  (14) 

V"  9  P 
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Figure  2.  Schematic  of  Computational  Cell  Structure 


FOR  U -NODES 


where  the  control  volume  around  the  grid  node  P  is  denoted  by  V*. 
In  Eq.  (14)  b  and  c  are  in  general  functions  of  ?.  The  represen¬ 
tation  of  the  convective  and  diffusive  fluxes  requires  some  elab¬ 
oration,  especially  since  present  research  has  implemented  some 
improvements  in  the  original  formulation  of  the  TEACH  code.  We 
discuss  this  aspect  in  detail  later  but  indicate  now  the  general 
form  which  represents  Eq.  (1),  after  the  flux  values  at  the 
control  volume  faces  are  approximated  in  terms  of  $  values  at 
the  neighboring  nodes.  For  the  control  volume  under  con¬ 
sideration  we  have 


(ap-b)  ,p  =  _  an;n  +  c,  (15) 

n 

where  ap  =  an,  L  denotes  the  summation  over  the  neighboring 
n  ‘  n 

nodes  N,  E,  S  and  W,  and  the  an ' s  and  ^  ’  s  are  respectively 
af.},...,  aw  and  .  .  .  ,  <p^.  Note  that  the  an '  s  represent  the 
coefficients  for  the  combined  convective  and  diffusive  fluxes 
of  4>. 

Equations  similar  to  Eq.  (15)  can  be  written  for  each  of  the 
dependent  variables  ?  at  every  node  P  in  the  interior  of  the  com¬ 
putational  domain.  Appropriate  modifications  to  the  expressions 
for  the  fluxes  are  necessary  when  the  nodes  adjoining  the  boun¬ 
daries  of  the  domain  are  encountered  (see  Paragraph  2.5).  The 
resulting  system  of  equations  can  be  solved  for  the  unknown  4)  by 
an  iterative  procedure,  if  the  needed  boundary  data  are 
available.  A  brief  discussion  of  the  boundary  conditions  is  pre¬ 
sented  in  Paragraph  2.5.  The  iterative  procedure  yields  an 
acceptable  solution  when  a  certain  convergence  criterion  is 
satisfied.  Although  this  question  is  considered  when  we  discuss 
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the  solution  algorithms  in  Paragraph  2.6,  we  point  out  that  a 
sufficient  condition  for  convergence  for  each  ?  is  that 


aD-b 


(  16a  ) 


for  all  equations,  and 


aP-b  !  >  l 

n 


(  15b) 


for  at  least  one  equation.  These  conditions  are  satisfied  if  all 
the  an ' s  are  positive  (since  l  an  =  ap)  and  b  is  neqative.  The 
latter  requires  that  the  linearization  in  Eq.  (14)  be  such  that  b 
is  negative. 

The  pressure  field  is  determined  from  an  equation  obtained 
by  combining  the  continuity  and  momentum  equations.  This  proce¬ 
dure  yields  a  finite-difference  pressure-correction  equation 
similar  to  Eq.  (15),  with  c  now  representing  the  local  mass 
imbalance  in  the  prevailing  velocity  field. 

We  now  return  to  the  question  of  evaluating  the  convective 
and  diffusive  fluxes  of  Following  is  a  discussion  of  dif¬ 

ferencing  schemes. 

The  derivation  of  the  fluxes  in  the  TEACH  code  can  be  better 
understood  by  considering  the  steady,  one-dimensional  equation 
for  the  convection-diffusion  balance  described  by 


d  ,  rT  .  d  ,  .  d  i, 

dfz  (  pwp)  ~  d7  (  r  d7} 


Applying  this  to  the  one-dimensional  grid  points  and  control 
volume,  we  obtain 


(  pW  ? )  -  (  pW  ? ) 

0  w 
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where  the  subscripts  e  and  w  respectively  denote  east  and  west 
faces  of  the  control  volume.  If  j  and  d j/dz  are  expressed  at 
control  volume  faces  in  terms  of  the  nodal  values  of  ?,  i.e.,  j^, 
jo  and  Jr,  Eq.  (13)  leads  to  the  algebraic  equation  for  the 
unknown  jp,  given  by  Eq.  (13). 

2.4.1  Exponential  Scheme 

A  good  starting  point  of  the  face-value  approximation  is  an 
analytical  solution  of  Eq.  (17),  which  is  easily  obtained  when 
jW  =  constant  and  r  =  constant.  Applying  the  analytical  solution 
to  the  intervals  E-P  and  W-P,  we  obtain  the  fluxes  into  the  right 
and  left  boundaries  of  the  control  volume,  given  respectively 

by, 33 

Je  ’  (pW>  "  rli}e  =  1  pW)e  UP+(  W7  ^xp(Pee)_1  -  (19a) 

and 

Jw  5  (pW*  -  rii)w  =  (?W)wi%+(  ^r?P)/>xp(Pew)_1  ^  '  (19b) 

where  Pee  and  Pew  are  the  local  Peclet  numbers.  These  are 
defined  as  Pee  =  (  pW )  e  (  -SzpE )  /  re  and  Pew  =  (  pW)w  (  5zwp )  /  rw  (  5z 
denotes  the  axial  grid  interval  between  adjacent  nodes).  Mote 
that  for  nonuniform  p  and  r,  average  values  must  be  used. 
Expressing  the  local  Peclet  numbers  in  terms  of  the  convective 
and  diffu  .ve  coefficients  C  5  pW  and  D  =  r/oz,  we  obtain 
Pee  =  Ce/De  an<3  Pew  =  Cw/Dw.  Through  the  use  of  Eqs.  (19)  and 
the  foregoing  definitions,  Eq.  (18)  can  be  cast  as 

ap  4>p  =  ar  ?£  +  Sir-i  ,  (  20  ) 


where  aE  h  Ce/  [exp( Pee ) -1 j ,  aw  5  Cwe  ;o l Pew ) /  ’exp ( Pew ) -1  ]  and 
ap  i  aE+aT/j.  In  obtaining  the  last  identity,  use  is  made  of  the 
equality  Ce=Cw  in  view  of  the  one-dimensicnal  mass  conservation. 
Note  that  Eq.  (20)  is  just  the  one-dimensional  analog  of  Eq. 
(15),  without  the  source  terms  b  jp-'-c . 


The  foregoing  development  describes  the  exponential  scheme 
for  the  one-dimensional  problem  where  it  provides  the  exact 
solution.  In  spite  of  its  highly  desirable  features,  however , 
the  exponential  scheme  has  not  been  widely  used  in  numerical 
computations,  because  the  exponentials  are  expensive  to  compute. 
Moreover,  the  scheme  dees  not  apply  exactly  in  two-  or  three- 
dimensional  situations,  and  for  nonzero  source  terms.  Therefore, 
the  extra  computational  expense  inherent  in  the  exponential 

scheme  does  not  justify  its  extension  to  the  two-dimensional 
problem  at  hand. 

2*4.2  "Hybrid"  Upwind  Scheme 

Accordingly,  the  developers  of  the  TEACH  code  have  exploited 
a  "hybrid"  differencing  scheme,  in  which  the  face  values  of  the 
fluxes  are  approximated  as  follows:39 

(Cw/2+V  V(V2-V  V  for  "2  <  Pew  <  2 
CwV  for  Pew  >  2 

CwV  for  Pew  <  -2 

(  21 ) 

'ce/2+De )  ^p  +  (ce/2-De )  ?E ,  for  -2  <  Pe0  <  2 
cefp«  for  Peo  >  2 

Ce  JE'  for  Pep  4  ~ 2  • 

Vie  note  that  in  obtaining  Eq.  (21),  the  exponentials  in  Eqs . 

(19)  have  been  approximated  by  piece-wise  linear  representations. 
The  resulting  expressions  for  the  fluxes  in  Eq.  ( 21)  correspond 
to  the  central  differencing  when  -2  <  Pe  <  2  and  the  upwind  dif¬ 
ferencing  (with  the  diffusive  contributions  being  neglected)  when 
?e  !  ,>  2,  thus  the  designation  "hybrid."  in  terms  of  the 


standard  formulation  of  Eq.  (20),  the  coefficients  for  the 
combined  convective  and  diffusive  fluxes  in  the  "hybrid"  upwind 
scheme  become 
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A  comparison  of  the  exponential  scheme  and  the  "hybrid" 
upwind  scheme  is  given  by  Patankar40  in  graphic  form.  This 
comparison  is  better  facilitated  by  the  introduction  of  the 
nond imens ional  coefficients  aE/0e  and  aw/Dw.  He  obtain 


aE/De  =  Peg/ vexp(Pee  )-l 

for  the  exponential  scheme  and 


l-Peo/2,  for  -2  <  Pep  <  2 

0  ,  for  Pe  >2 

e 

-Pe  ,  for  Pe  *  -2 

s  0 


(  22a) 


(  22b) 


for  the  "hybrid"  upwind  scheme 


The  use  of  the  ’’hybrid"  upwind  differencing  scheme  has  been 
generally  accepted  by  the  TEACH-code  modelers  and  little  loss  of 
accuracy  has  been  claimed  in  following  this  scheme.  It  should  he 
noted,  however,  that  the  "hybrid"  scheme  yields,  according  to  Eg. 
(22b),  a  value  of  0  at  ?ee  =  2  and  a  value  of  2  at  ?e0  =  -2  for 
ar/De.  The  exact  values  from  the  exponential  scheme,  according 
to  Eq.  (22a),  are  respectively  0.3130  and  2.313.  Thus,  at 
|  Pee  j  =  2  the  departure  of  the  "hybrid"  upwind  solution  from  the 
exact  (exponential)  solution  of  the  one-dimensional  convection- 
diffusion  equation  is  rather  large. Furthermore,  it  is  clearly 
arbitrary  to  neglect  the  diffusion  effects  as  soon  as  !  ?e  | 
exceeds  2  _ ( see  Eq.  (21)  . 

2.4.3  Power-Law  Differencing  Scheme 

A  better  approximation  to  the  exact  solution  has  been 
obtained  by  the  power-law  differencing  scheme,  which  is  described 
in  Patankar.4!  A  comparison  of  the  "hybrid"  upwind  and  power-law 
schemes  in  Reference  40  shows  that  the  solutions  become  identical 
for  |  Pe  |  >  10.  For  example,  the  nondimens ional  coefficient 
aS/De  f°r  ?ower-law  scheme  is  given  by 

-Pe  ,  for  Pe  <  -10 

e  e 

(1+0.1  Pe  ) 3 -Pe ,  for  -10  <  Pe  <0 
e  e 

,  (22c) 

(1-0.1  Pe  ) 3 ,  for  0  <  Pe  <10 
e  e 

0  ,  for  Pe  >  10  . 

e 

From  Eq.  (22c)  it  is  obvious  that  the  power-law  scheme  solutions 
are  much  closer  to  the  exponential-scheme  solutions  than  the 
"hybrid"  upwind  solutions,  for  -10  <  Pee  <  0  and  for  0  < 

Pee  <  10.  We  noted  earlier  that  at  j  Pee  |  =  2,  the  "hybrid" 
upwind  solution  is  the  farthest  from  the  exponential  solution. 

On  the  other  hand,  Eq.  (22c)  yields  that  a;r/De  is  2  .  328  at 
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?ee  =  -2  and  is  T .  3  2”'"’  a  t  °ee  =  2.  Thus,  the  power-law  scheme 
provides  an  extremely  zoo<i  representation  of  the  exponential 
oehavior  f or  t re  3e:le'. — .  umcer  range  where  the  "hybrid"  scheme  is 
not  s a  t  :  a f  a  c t '  r v . 

‘rm  tne  *  're  :o  :n  a  comparisons,  it  appears  that  for  the 
one— 1  men  s  :  nn  a  1  iiif  us:  on -convection  equation  the  power-law 
scheme  :s  superior  to  the  "hybrid"  upwind  scheme.  Moreover,  the 
diffusive  effects  are  retained  in  the  power-law  scheme  for  a 
larger  Recler-number  range,  viz.,  |  Pe  }  <•  10.  In  view  of  these 
advantages  and  since  tne  power-law  expressions  are  not  par¬ 
ticularly  expensive  to  compute,  we  have  employed  the  two- 
dimensional  extension  of  the  power-law  scheme  {due  to 
Patan’xar^'  in  our  calculations  with  the  TEACH  code.  As  part  of 
our  numerical  computations  of  the  centerbodv  combustor 
flowfields,  we  have  compared  the  predictions  of  the  "hybrid" 
upwind  and  power-law  schemes.  This  aspect  is  discussed  in 
Section  4. 

2.5  BOUNDARY  CONDITIONS 

The  solution  of  Eq.  (1)  requires  the  specification  of 
appropriate  boundary  conditions  for  each  of  the  dependent 
variables.  Since  the  governing  equation  is  elliptic,  the 
boundary  conditions  must  be  prescribed  on  all  the  boundaries  of 
the  computational  domain  shown  schematically  in  Figure  1.  Note 
that  because  of  symmetry  only  the  top  half  of  the  centerbody 
combustor  configuration  is  represented  by  the  computational 
domain . 

The  top  boundary  of  the  computational  domain  is  the  duct 
wall  and  the  bottom  boundary  is  the  axis  of  symmetry.  The  left 
boundary  represents  the  inflow  boundary  and  consists  of  the 
annular  air  inlet,  the  centerbody  face,  c.  .  the  central  fuel 
inlet.  The  right  boundary  is  the  outflow  boundary,  the  location 
of  which  is  unknown  a  priori.  Therefore,  the  specification  of 
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the  location  of  this  boundary  is  essentially  arbitrary  and,  as 
discussed  in  Section  1,  it  is  essential  to  ensure  that  the 
sensitivity  of  the  computed  solutions  to  the  specification  of 
this  boundary  location  is  not  significant. 

The  boundary  conditions  selected  in  our  computations  are 
shown  in  Table  2.  Some  general  observations  can  be  made  with 
respect  to  the  conditions  at  the  four  boundaries.  Along  the 
bottom  boundary  representing  the  centerline,  symmetry  con¬ 
siderations  require  that  3  $/  3r  =  0  for  all  "j>;  in  addition, 
because  of  symmetry,  the  radial  velocity  U  vanishes  on  the  cen¬ 
terline.  Along  the  top  boundary  representing  a  rigid  impermeable 
wall,  the  velocity  components  are  set  equal  to  zero  (due  to  the 
no-slip  requirement  on  W  and  the  nonmoving  wail  for  rJ)  .  However, 
to  minimize  the  computer  storage  and  run  time,  the  dependent 
variables  at  the  wall  are  linked  to  those  at  the  grid  node 
adjacent  to  the  wall  by  equations  which  are  consistent  with  the 
logarithmic  law  of  the  wall.  For  the  tangential  velocity  (W  in 
this  case),  this  results  in  a  condition  on  the  wall  shear  stress 
tw  in  the  calculations.  Wall-function  formulations  also  govern 
the  boundary  conditions  of  k  and  e.  For  the  species  mass 
fraction  Yp,  the  impermeable  boundary  requires  the  normal 
derivative  (  3Yp/ 3r  at  the  top  boundary)  to  vanish.  In  the  code's 
finite-difference  formulation,  this  is  accomplished  through 
setting  the  coefficient  a^  to  zero  for  the  cell  whose  north  boun¬ 
dary  coincides  with  the  top  boundary.  Along  the  left  boundary, 
similar  considerations  apply  to  the  rigid  impermeable  boundary 
representing  the  centerbody  face:  for  the  U  velocity  the  law  of 
the  wall  is  employed  to  express  tw.  For  k  and  £  the  wall  func¬ 
tions  are  used.  The  vanishing  of  3Yp/ 3z  along  this  boundary  is 
implied  by  setting  the  coefficient  a^  to  zero  for  the  cells  with 
the  west  boundaries  coinciding  with  the  centerbody  face.  Along 
the  inlet  ports  in  the  left  boundary,  the  conditions  are  spe¬ 
cified  through  experimental  measurements.  The  inlet  profiles  of 
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TABLE  2 

BOUNDARY  CONDITIONS 
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W  are  obtained  from  experiments.  For  well-designed  inlets,  the 
radial  velocity  U  is  typically  close  to  zero.  Any  nonzero  radial 
velocity  distribution,  if  specified,  must  satisfy  the  continuity 
equation  at  the  inlet.  The  inlet  distributions  of  species 
mass  fraction  Yp  are  usually  assumed  to  be  uniform  across  the 
respective  (central)  fuel  and  (annular)  oxidant  ports. 

The  distributions  of  k  and  e  remain  to  be  specified  at  the 
inlet.  The  specification  of  k  entails  the  prescription  of 
appropriate  value  of  the  parameter  TURBIN  (FORTRAN  variable  in 
the  TEACH  code).  The  inlet  profile  of  k  is  obtained  from 

k  =  TURBIN  x  W2  ,  (23) 

m 

where  Wj_n  denotes  the  mean  axial  velocity  at  the  inlet.  A  value 
of  0.03  is  considered  typical  for  TURBIN.  Note  that  experimental 
data  on  the  turbulence  intensity  measurements  in  the  inlet  pro¬ 
vide  a  basis  for  selecting  the  appropriate  value  of  TURBIN.  Of 
course,  in  the  absence  of  the  turbulence  intensity  measurements 
in  all  three  orthogonal  coordinate  directions,  recourse  must  be 
made  to  the  assumption  of  isotropy  in  determining  k  from  the 
turbulence  intensity  results  in  one  or  two  directions. 

The  specification  of  the  inlet  distributions  of  e  is  not  yet 
possible  from  the  measurement  of  the  dissipation  rate.  As 
discussed  in  Paragraph  2.2,  it  becomes  necessary  to  introduce 
certain  assumptions  on  the  turbulence  length  scale  in  the  inlet. 
Depending  on  the  choice  of  the  length  scale  Z  ^  [see  Eq .  (7b)’ 
or  1 2  see  Eq .  (9b)]  ,  appropriate  values  of  the  parameter  X 
(denoted  bv  the  FORTRAN  variable  ALAMDA  in  the  TEACH  code)  must 
be  selected.  Table  3  shows  some  examples  of  the  length  scales 
used  by  different  authors.  It  is  clear  that  the  length  scales  in 
different  flow  configurations  differ  considerably.  This  obser¬ 
vation  should  hardly  be  surprising,  since  no  single  length  scale 
of  turbulence  can  be  expected  to  represent  all  types  of  turbulent 
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TABLE  3 

SPECIFICATION  OF  c  AT  THE  INLET 
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flows.  What  is  crucially  important,  however,  is  to  determine  how 
strongly  the  flowfield  predictions  depend  on  the  inlet  length 
scale.  Our  numerical  results  discussed  in  Sections  III  and  IV 
address  this  question. 

For  the  centerbody  configuration,  the  rate  of  dissipation  at 
the  inlet  is  given  by  Eq.  (7b),  with  the  length  scale  given  by 

1^=a5,  (24) 

where  6  is  a  characteristic  reference  length  (see  Table  3).  For 
the  central  port,  the  reference  length  is  the  tube  radius  Rg  . 

For  the  annular  port,  the  reference  length  is  the  so-called 
hydraulic  mean  radius  which  is  given  by  R2~R]_/  where  s2 
duct  radius  and  Ri  is  the  centerbody  radius.  We  note  that  the 
value  of  0.03  for  X  used  by  Sturgess  and  Sved*’  corresponds  to  a 
value  of  0.3333  when  Eqs .  (7a)  and  (7b)  are  used  (as  in  the 

present  study)  instead  of  Eqs.  (9a)  and  (9b). 

2.6  SOLUTION  PROCEDURE 

Writing  Eq.  (15)  for  each  node  in  the  domain  of  our 
interest,  systems  of  algebraic  equations  for  each  dependent 
variable  4)  are  obtained.  The  resulting  large  nonlinear  coupled 
systems  naturally  require  an  iterative  solution  procedure.  This 
involves : 

(a)  solving  the  system  of  equations  for  one  dependent 
variable  with  linearization  (inner  iteration),  and 

(b)  connecting  properly  the  inner  iteration  for  all  s 
until  the  converged  solution  is  obtained  (outer 

i teration )  . 
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The  outer  iteration  loop  consists  of  the  following 
sequence :  ^ 

(a)  guess  the  pressure  field  d, 

(b)  solve  the  momentum  equations  to  obtain  VI*,  u*, 

(c)  solve  the  pressure  correction  equation  (based  on  mass 

conservation)  to  obtain  the  correction  to  the  pressure 
field, 

(d)  update  W  and  U  based  on  the  pressure  correction  and 

★ 

update  p  to  get  W,  U,  and  p, 

(e)  solve  the  systems  for  other  $'s, 

(f)  treat  $  as  p*  and  return  to  step  (b),  and  repeat  the 

whole  procedure  until  convergence  is  achieved. 

The  inner  iteration  loop  employs  a  line-by-line  solution 
procedure,  for  which  the  well  known  Tri-digonal  Matrix  Algorithm 
(Thomas  Algorithm)  is  used.  For  each  <p,  the  number  of  sweeps  is 
specified.  One  sweep  is  composed  of  an  east-west  sweep  and  a 
north-south  sweep.  The  number  of  sweeps  for  each  variable  ?  is 
set  in  the  main  program  by  the  user. 

Currently,  the  number  of  sweeps  for  pressure  correction  is 
set  at  4,  while  that  of  all  the  other  variables  is  set  at  2. 

It  should  be  noted  that  the  specification  of  large  number  of 
sweeps  is  not  desirable,  since  the  solution  is  only  tentative 
after  all. 

2.6.1  Underrelaxation 

Underrelaxation  is  usually  employed  for  nonlinear  problems 
to  avoid  divergence  in  the  iterative  solution  process. 

For  all  the  dependent  variables  other  than  pressure,  the 
underrelaxation  comes  through  the  modification  of  the  coef¬ 
ficients  of  Eq.  (15),  that  is,  the  underrelaxed  equations  are 
solved  by  the  line-by-line  solver. 
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Let  a  be  the  underrelaxation  factor.  It  can  be  readily 
shown  that  the  underrelaxed  form  of  Eq.  (15)  is  given  as 


a 


a  j> 

n  n 


+  c 


(1-a) 


* 


where  ap  =  ap-b.  This  is  the  equation  which  is  actually  solved 
by  the  line-by-line  solver. 

The  pressure  density  and  viscosity  (effective  viscosity)  are 
underrelaxed  via 


?r.ew  ~  afnew  +  ^  ^  ?old‘ 

In  the  present  calculations,  all  the  underrelaxation  factors  are 
set  at  0.5. 

2.6.2  Convergence 

The  convergence  criterion  is  somewhat  arbitrary  and  is 
supplied  by  the  user.  The  present  code  uses  the  discretized 
Eq.  (15)  to  derive  the  convergence  criterion.  From  Eq.  (15)  a 
residual  R  is  defined  as 

R  =  ^  an*n  +  c  '  (aP-b)  V 

n 

Hence,  if  the  discretized  equation  is  satisfied,  R  will  be  zero. 
Presently,  the  sum  of  all  )  R  |  *s  over  the  interior  grid  points  is 
normalized  by  appropriate  reference  value  for  each  ?,  and  the 
resulting  values  are  used  to  test  convergence.  All  the  residuals 
are  monitored  to  see  the  convergence  history.  Numerical  calcula¬ 
tions  reported  here  are  based  on  the  criterion  that  the  maximum 
of  the  normalized  residuals  be  less  than  10_<*  for  flows  involving 
a  single  species,  and  be  less  than  10“2  for  flows  involving  two 


spec ies . 


SECTION  3 

PRELIMINARY  NUMERICAL  COMPUTATIONS 


As  indicated  in  Paragraph  1.3,  the  numerical  modeling  of  the 
centerbodv  combustor  flowfields  represented  several  phases  and 
addressed  a  number  of  aspects.  These  ranged  from  preliminary 
computations  directed  toward  the  familiarization  and  examination 
of  the  AFWAL/PORT  version  of  the  TEACH-T  program  to  full-fledged 
numerical  simulation  and  comparison  with  experimental  data.  A 
summary  of  the  preliminary  calculations  completed  in  the  oresent 
modeling  efforts  is  presented  in  Table  4.  In  this  section  we 
discuss  the  details  of  these  preliminary  computations.  Further 
modeling  refinements  and  selected  results  are  discussed  in  the 
next  section. 

3.1  APL  CONFIGURATION 

The  initial  computations  for  the  isothermal  flowfield 
with  2  kg/s  air  and  8  kg/hr  CO2  flows  resulted  in  numerical 
instabilities,  convergence  failures,  and  nonphysical  flowfields. 
These  computations  used  the  two-equation  turbulence  model  in 
the  code.  Careful  examination  of  the  code  revealed  that  the 
AFWAL/PORT  version  of  the  code  was  written  only  for  the  reactino- 
flowfield  computations.  At  the  suggestion  of  Dr.  Roquemore,  we 
contacted  Mr.  Russ  Claus^  0f  NASA/Lewis  who  had  been  operating  a 
version  of  the  "TEACH"  Code  to  predict  the  APL  combustor 
flowfields.  It  turned  out  that  this  version  was  specifically 
written  for  nonreacting-f lowf ield  computations  involving 
one-component  fluid. 

Accordingly,  we  decided  to  modify  the  AFWAL/PORT  version  to 
perform  the  nonreacting-f lowf ield  computations.  For  this 
purpose,  we  resorted  to  the  artifice  of  specifying  zero  heat  of 
formation  for  the  "fuel"  (i.e.,  COg).  This  resulted  in  the 
desired  isothermal  flowfield  as  evidenced  by  the  computed 
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temperature  field.  Also,  the  value  of  the  underrelaxation 
factors  (see  Paragraph  2.6.1)  for  all  the  dependent  variables  was 
changed  to  0.5  to  ensure  better  convergence.  However,  the 
computed  fields  of  axial  velocity  and  CO2  mass  fraction  revealed 
nonphysical  characteristics.  Further  examination  of  the  code 
became  necessary  to  determine  the  source  of  these  uncertainties. 
We  found  that  the  AFWAL/PORT  version  requires  correct 
specification  of  excess  air  to  determine  the  air  flow  rate  from 
the  specified  "fuel"  flow  rate.  Thus,  even  when  nonreacting 
calculations  are  performed,  depending  on  the  flow  rates  of  air 
and  CO2  used  in  the  computations,  the  appropriate  value  of  the 
excess  air  must  be  determined  and  furnished  as  input  to  the 
program.  Once  this  requirement  was  satisfied,  the  mass-flow 
calculations  became  consistent  with  the  inlet  velocity  profiles 
and  the  numerical  computations  proceeded  smoothly. 

3.1.1  Variations  of  Annular  and  Central  Jet  Flow  Rates 

The  appropriately  modified  version  of  the  TEACH  Code  was 
used  to  compute  the  isothermal  flowfields  of  the  centerbody 
combustor  for  different  combinations  of  annular  and  central-jet 
flow  rates.  Thus,  the  completed  predictive  modeling  calculations 
correspond  to  (a)  zero  annular  flow  and  4,  8,  12,  and  16  kg/hr 
central  CO2  flow;  (b)  0.05  kg/s  annular  air  flow  and  4,  8,  12, 
and  16  kg/hr  central  CO2  flow;  (c)  0.5  kg/s  annular  air  flow  and 
0  (actually  10— 3 ) ,  1,  2,  4,  8,  12,  and  16  kg/hr  central  CO2  flow; 
(d)  1  kg/s  annular  air  flow  and  0,  2,  4,  6,  8,  12,  and  16  kg/hr 
central  CO2  flow;  (e)  2  kg/s  annular  air  flow  and  0,  4,  8,  12, 

16,  18,  20,  and  22  kg/hr  central  CO2  flow;  and  (f)  3  kg/s  annular 
air  flow  and  4,  6,  8,  12,  16,  and  18  kg/hr  central  CO2  flow. 

In  (b)  above,  calculations  were  also  made  for  1  x  10”^  kg/s 
annular  air  flow  and  4  kg/hr  CO2  flow.  These  calculations 
yielded  results  identical  with  those  of  the  zero  annular  flow  and 
4  kg/hr  central  flow.  Finally,  reacting  flowfield  modeling  was 
done  for  2  kg/s  air  flow  and  propane  flow  in  the  central  jet. 
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Only  the  case  of  2  kg/hr  propane  resalted  in  converged  3'lations 
that  were  physically  acceptable. 

All  these  modeling  calculations  were  conducted  for  the 
41-X-34  computational  grid  with  nonuniform  grid  spacing  in  both 
axial  and  radial  directions.  The  axial  extent  of  the 
computational  domain  was  0.3m.  In  both  the  annular  and 
central-jet  inlets,  only  uniform  axial  velocity  profiles  ’were 
used.  In  the  two-equation  turbulence  model,  standard  values  for 
the  various  "constants"  were  used.  Finally,  for  the 
specification  of  the  values  of  the  turbulent  kinetic  energy  and 
the  dissipation  rate  (through  the  specification  of  the  inlet 
length  scale),  only  the  default  values  (TURBIN  =  0.03  and  ALAMDA 
=  0.005)  in  the  code  were  used. 

The  computed  flowfields  for  the  two-jet  flows  (2  kg/s 
annular  air  flow  and  several  CO2  central-jet  flows)  exhibited 
trends  that  were  in  conformity  with  the  earlier  Laser  Doppler 
Anemometer  measurements.  Quantitatively,  however,  the 
calculations  underpredicted  the  centerline  forward 
stagnation-point  distances  and  overpredicted  the  rear 
stagnation-point  distances.  It  appeared  that  better  agreement 
between  the  predicted  and  measured  results  could  be  obtained  by 
varying  the  inlet  turbulence  length  scales. 

Some  of  the  calculated  results  were  plotted.  These  included 
the  axial  profiles  of  the  centerline  mean  velocities,  turbulence 
intensities,  and  CO2  mass  fractions.  The  radial  profiles  of  the 
axial  velocity  at  an  axial  location  of  12  cm  from  the  centerbody 
face  were  plotted  for  the  central  flow  rates  of  0,  8,  and 
16  kg/hr.  All  these  profiles  showed  the  trends  observed 
experimentally.  The  axial  profiles  were,  with  respect  to  the 
axial  distance,  normalized  by  the  centerbody  diameter.  The 
centerline  mean  velocity  profiles  were  also  obtained  by 
normalizing  the  axial  distance  with  the  centerline  rear 
stagnation  point  distance.  However,  the  quantitative  agreement 
between  the  prediction  and  measurement  was  poor. 
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The  single-jet  calculations  (with  zero  annular  flow)  showed 
that  the  results  correspond  to  those  of  a  free  jet.  The  cen¬ 
terline  axial  velocity  profiles  (velocity  normalized  by  the  jet 
exit  velocity  and  the  axial  distance  normalized  by  the  initial 
jet  diameter)  exhibited  on  a  log-log  plot  clearly  indicated  the 
free  turbulent-jet  behavior.  The  results  for  the  four  central 
flow  rates  of  4,  8,  12,  and  16  kg/hr  correspond  to  a  single  curve 
which  showed  a  potential  core  up  to  5  diameters,  a  slope  of  -1 
between  7  and  55  diameters,  and  an  axial  velocity  decaying  to  10 
percent  of  the  initial  velocity  in  55  diameters.  The  axial 
variation  of  the  centerline  turbulent  intensity  was  also  plotted. 
When  the  turbulence  intensity  was  normalized  with  respect  to  the 
initial  jet  velocity,  the  axial  profile  showed  an  initial  decay 
in  the  potential  core,  a  rapid  rise  between  5  and  10  diameters, 
and  a  rapid  decay  beyond.  However,  when  the  normalization  was 
with  respect  to  the  local  values  of  the  centerline  mean  axial 
velocity,  there  was  a  slow  decay  (from  3  to  10  percent  to  2  to  4 
percent)  up  to  four  diameters,  a  rapid  rise  to  30  percent 
(between  10  and  20  diameters),  and  a  gradual  rise  to  32  percent 
beyond.  While  it  is  not  clear  if  the  initial  decay  is  real,  the 
rapid  rise  and  the  plateauing  are  consistent  with  the  earlier  LDA 
results.  Because  of  the  boundary  condition  representing  the  con¬ 
fined  duct,  in  the  absence  of  annular  air  flow  the  centerline 
CO2  mass  fraction  did  not  show  any  significant  decay  downstream 
of  the  jet  exit.  Calculations  for  4  kg/hr  CO2  flow  with  1  x 
10~5  fcg/s  annular  air  flow  showed  identical  results. 

To  examine  the  concentration  decay  at  small  annular  air 
flows,  two-jet  calculations  were  made  for  0.05  kg/s  annular  air 
flow  and  4,  8,  12,  and  16  kg/hr  CO2  flow.  The  computed  results 
for  this  set  also  showed  the  free-jet  behavior  for  the  centerline 
axial  velocity  (normalized  with  respect  to  the  initial  velocity) 
and  the  centerline  CO2  mass  fraction.  Indeed,  both  profiles 
appeared  to  be  identical  (with  the  results  of  the  four  flow  rates 
coinciding ) . 
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The  reacting  flowfield  calculat  ions  for  2  <3 /hr  propane  flo 
showed  that  the  centerline  variation  of  near,  axial  velocity,  tar 
bulent  intensity,  CO 2  concentration ,  and  temperature  was  con¬ 
sistent  wtih  the  expected  trends.  Hinder  central  flow  rates  din 
not  yield  convergent  results. 

Our  preliminary  experience  with  the  TEACH  Code  indicated 
reasonably  correct  trends  in  predicting  the  centerbodv  combustor 
flowfields.  Further  computations  of  the  isothermal  flov/fields 
with  different  grid  spaci.nqs,  with  nonuniform  inlet  velocity  pro 
files,  and  with  appropriate  inlet  length  scales  were  necessary 
before  improvements  in  quantitative  predictions  and  scaling  cri¬ 
teria  could  be  obtained.  Accordingly,  questions  relating  to  the 
grid  i.odeoendence  of  the  computed  solutions,  boundary-layer 
deve  looment  in  the  inlets,  and  inlet  turbulence  in.  ten  si  tv  and 
length  scales  were  addressed  in  several  sensitivity  tests. 

3.1.2  Sensitivity  Tests 

Our  calculations  with  2  kg/s  air  flow  rate  and  4,  3,  and 
12  kg/hr  CO2  flow  rates  examined  each  of  the  above  questions 
separately,  as  discussed  in  the  following  paragraphs. 

3.1.2.  1  Grid  Independence  of  Solutions 

Numerical  experiments  were  carried  out  with  the  uniform 
inlet  velocity  profiles  and  default  values  for  the  inlet  kinetic 
energy  and  length  scale  parameters;  however,  nonuniform  grid 
intervals  50  percent  larger  than  the  earlier  ones  were  employed. 
The  results  of  the  centerline  stagnation  points  (distances  in 
meters)  are  shown  in  Table  5. 

The  earlier  results  with  finer  grid  spacmgs  had  underpre¬ 
dicted  the  locations  of  the  forward  stagnation  points  and 
overpred icted  those  of  the  rear  stagnation  ooints.  The  use  of 
cruder  grid  intervals,  however,  decreases  both  stagnation  dis¬ 
tances,  thereby  making  the  disagreement  between  measurement  and 
prediction  for  the  forward  stagnation  distance  worse  (except  the 
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TABLE  5 

INFLUENCE  OF  GRID  SPACINGS 
ON  (APL)  CENTERLINE  STAGNATION  DISTANCES 


% 


CENTRAL  JET 
FLOW  RATE 
Kg/hr 


I  ENTERLINE  STAGNATION  DISTANCES.  Ti 


FORWARD 

STAGNATION 

REAR  STAGNATION 

1  INITIAL  GRID 

50%  LARGER 
(CRUDER) 

INITIAL  GRID 

50%  LARGER 

0.0198 

0.0185 

0.184 

0.176 

0.0562 

0.0548 

0.139 

0.130 

12 


0.0854 


0.0871 


0.191 


0.133 


12  kg/hr  case).  All  the  same,  the  changes  due  to  the  cruder  grid 
intervals  are  only  on  the  order  of  five  percent.  Thus,  it  would 
seem  that  the  computations  do  exhibit  adequate  grid  independence 
in  the  range  of  grid  intervals  investigated  in  the  computations. 

3. 1.2. 2  Influence  of  Inlet  Velocity  Profile 

The  effect  of  the  boundary-layer  development  in  the  inlets 
is  shown  in  Table  6.  For  this  purpose,  the  inlet  velocity  pro¬ 
files  employed  were  uniform  everywhere  except  at  the  grid  nodes 
adjacent  to  the  inlet  walls  where  the  velocities  were  95  percent 
of  the  uniform  values.  The  older  grid  intervals  and  default 
values  of  TURBIN  and  ALAMDA  were  used. 

The  boundary-layer  growth  in  the  inlets  tends  to  increase 
both  stagnation  distances.  VJhile  this  is  in  the  desired  direc¬ 
tion  for  the  forward  stagnation  distance,  it  makes  the  departure 
from  the  measurement  for  rear  stagnation  distance  worse. 

However,  the  change  in  the  latter  case  is  only  on  the  order  of 
three  percent,  thus  showing  that  the  effect  on  the  rear  stagna¬ 
tion  point  is  very  small.  In  the  case  of  the  forward  stagnation 
point,  the  change  is  about  seven  percent  for  4  kg/hr,  four  per¬ 
cent  for  3  kg/hr,  and  five  percent  for  12  kg/hr.  It  may  be 
recalled  that  the  departure  of  the  prediction  from  the  measure¬ 
ment  is  also  higher  at  lower  central  flow  rates.  Thus,  it 
appears  that  the  use  of  the  nonuniform  inlet  velocity  profiles 
would  result  in  better  agreement  between  measured  and  predicted 
forward  stagnation  distances,  without  causing  significantly 
greater  disagreement  for  the  rear  stagnation  distances. 

3. 1.2. 3  Influence  of  Inlet  Turbulence  Parameters 

Computations  that  examined  the  influence  of  inlet  turbulent 
kinetic  energies  and  length  scales  employed  the  values  of  0.06 
for  TURBIN  and  0.05  for  ALAMDA  ( X).  These  are  respectively  twice 
and  ten  times  the  default  values  for  these  parameters.  Flat 
inlet  velocity  profiles  and  finer  grid  intervals  were  used. 

Table  7  presents  the  results. 
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TABLE  6 

INFLUENCE  OF  INLET  VELOCITY  PROFILES 
ON  (APL)  CENTERLINE  STAGNATION  DISTANCES 


! 

I  CENTRAL  JET  I _ CENTERLINE  STAGNATION  DISTANCES,  m 


FLOW  RATE 

FORWARD 

STAGNATION 

REAR  STAGNATION 

)^g/hr 

FLAT  PROFILE 

FLAT  PROFILE 

95%  PROFILE 

ii i 

4 

| 

0.0198 

0.0212 

0.184 

0.189 

¥ 

8 

0.0562 

0.0582 

0.189 

0.193 

2 


o.oas-t 


0.0898 


0.191 


0.196 
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It  is  clear  that  the  effect  on  the  rear  stagnation  distance 
is  on  the  order  of  only  one  percent.  For  the  forward  stagnation 
distance,  the  change  decreases  from  about  3.5  percent  at  the 
lowest  central  flow  rate  to  about  0.8  percent  at  the  highest. 
Again,  the  perceived  chanqes  favor  the  forward  stagnation  dis¬ 
tances.  To  isolate  the  influence  of  the  length  scales  from  that 
of  the  turbulent  energies,  calculations  were  made  for  different 
values  of  TURBIN  and  ALAMDA  for  the  central  flow  rate  of  3  kg/hr. 
These  results  are  shown  in  Table  3. 

It  is  seen  that  doubling  the  length  scale  increases  the  for¬ 
ward  stagnation  distance  by  0.9  percent  and  the  rear  stagnation 
distance  by  0.5  percent.  However,  doubling  the  turbulent  kinetic 
energy  (although  at  the  larger  length  scales)  increases  the  for¬ 
ward  stagnation  distance  by  1.4  percent  and  the  rear  stagnation 
distance  hy  1.6  percent.  Doubling  the  kinetic  energy  and 
increasing  the  length  scales  by  ten  times  appear  to  be  better 
tljan  doubling  both  the  kinetic  energy  and  the  length  scales  in 
the  sense  that,  in  the  former  case,  the  increase  in  forward 
stagnation  distance  is  accompanied  by  a  lesser  increase  in  rear 
stagnation  distance.  However,  it  must  be  reiterated  that  the 
perceived  changes  are  too  small  to  be  significant,  especially  at 
higher  central  flow  rates. 

3. 1.2. 4  Influence  of  the  Exit  Boundary  Location 

A  single  computation  for  the  central  flow  rate  of  4  kg/hr 
was  made  with  flat  inlet  velocity  profiles  and  default  values  of 
ALAMDA  and  TURBIN  but  with  an  axial  computational  extent  of 
1.2  m.  This  corresponded  to  the  39-X-39  computational  grid  of 
Sturgess  and  Syed . ®  Since  the  radial  extent  is  fixed  by  the  duct 
radius  (of  12.7  cm),  the  radial  grid  intervals  were  kept  the  same 
as  our  earlier  (finer)  grid  but  the  axial  grid  points  were  those 
of  the  39-X-39  grid.  This  resulted  in  a  forward  stagnation 
distance  of  0.0209  m  and  a  rear  stagnation  distance  of  0.13  m  (as 
compared  to  our  benchmark  values  of  0.0198  m  and  0.134  m).  To 
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TABLE  8 

INFLUENCE  OF  ALAMO  A  .AND  TJRBXN  FOR  3  <g/hr  COj  FLOW 


TURBIN 
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"iNTERLINE  I 
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0.03 

0.285 
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0.139 

0.03 

0.57 

0.0567 
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0.06 

0.57 
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investigate  this  computational  grid  further,  flowfield  com¬ 
putations  were  completed  for  the  three  central  jet  CO2  flow  rates 
of  4,  3,  and  12  kg/hr.  Nonuniform  inlet  velocity  profiles  of  the 
type  used  earlier  and  the  values  of  0.05  and  0.06  respectively 
for  ALAMDA  and  TURBIN  were  employed  in  these  calculations.  The 
centerline  distances  of  the  forward  and  rear  stagnation  points 
are  as  follows:  4  kg/hr  :  0.0216  m  and  0.188  m;  3  kg/hr  : 

0.0634  m  and  0.195  m;  12  kg/hr  :  0.0964  m  and  0.196  m.  Thus,  it 
became  clear  that  all  these  computations  predicted  centerline 
rear  stagnation  distances  typically  50  percent  larger  than  the 
measured  distances.-^  Perceived  influences  of  changes  in  grid 
intervals,  inlet  velocity  profiles,  inlet  turbulent  energies  and 
inlet  length  scales  did  not  result  in  better  agreement  between 
the  predicted  and  measured  values. 

3.2  UCI  CONFIGURATION 

A  centerbody  combustor  of  roughly  1/5  scale  model  of  the  APL 
configuration  is  being  operated  at  the  University  of  California, 
Irvine.  Since  some  experimental  data  were  available  for  this 
combustor, 21  it  was  of  interest  to  obtain  the  predictions  with 
our  TEACH  code.  For  this  purpose,  the  computations  have  employed 
values  of  50.8  mm,  28  mm,  and  1.27  mm  respectively  for  the 
diameters  of  the  duct,  the  centerbody,  and  the  fuel  tube. 
Subsequently,  we  were  told  that  the  actual  UCI  dimensions  were 
slightly  different,  viz.,  51  mm,  30.5  mm,  and  1.3  mm.  It  is 
noted  that  neither  set  exactly  corresponds  to  a  1/5  scale  model 
of  the  APL  combustor.  Our  computations  were  carried  out  with 
respect  to  the  former  set  of  dimensions. 

3.2.1  Variations  of  Annular  and  Central  Jet  Flow  Rates 

The  slight  difference  in  the  dimensions  of  the  configuration 
indicated  earlier  meant  that  a  reference  duct  velocity  of  7.5  m/s 
corresponds  to  an  annulus  velocity  of  11.7  m/s  in  the  UCI  experi¬ 
ments  and  an  annulus  velocity  of  10.77  m/s  in  our  calculations. 
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All  the  calculations  for  sensitivity  tests  were  done  for  the 
annulus  air  velocity  of  10.77  m/s.  The  UCI  experimental  results 
were  reported  in  terms  of  the  overall  equivalence  ratios.  Thus, 
for  an  equivalence  ratio  of  0.02  and  annular  air  velocity  of 

10.77  m/s,  the  central  CO2  inlet  velocity  was  11.04  m/s.  For 
this  value  of  CO2  velocity  and  an  equivalence  ratio  of  0.05,  the 
annular  air  velocity  was  4.31  m/s.  Accordingly,  the  predictive 
calculations  corresponded  to  (a)  annular  air  velocity  of 

10.77  m/s  and  central  CO2  velocities  of  0,  5,  10,  15,  20,  25  and 
30  m/s;  (b)  annular  air  velocity  of  4.31  m/s  and  central  CO2 
velocities  of  0,  2.5,  5,  7.5,  10  and  15  m/s.  In  both  of  the 
above  cases,  the  highest  CO2  velocity  cases  resulted  in  the 
complete  penetration  of  the  recirculation  region  by  the  central 
jet. 

3.2.2  Sensitivity  Tests 

The  numerical  calculations  that  examined  the  sensitivity  of 
model  parameters  were  based  upon  the  annular  air  velocity  of 

10.77  m/s  and  central  CO2  velocity  of  11.04  m/s.  The  two  major 
tests  related  to  the  turbulence  Schmidt  number  a£  for  the  dissi¬ 
pation  equation  (see  Table  1)  and  the  constant  cp  in  the  k- t 
model  [see  Eq .  ( 4 )  ] . 

3.2.2. 1  Influence  of  cr£ 

In  the  TEACH  code,  the  eddy  viscosity  (the  momentum  exchange 
coefficient)  is  computed  from  Eq.  (4).  The  exchange  coefficients 
for  other  variables  are  obtained  by  prescribing  the  values  of  the 
appropriate  Schmidt  numbers,  a.  Thus,  as  shown  in  Table  1, 
constant  values  of  1  and  1.3  are  used  for  and  a£, 
corresponding  to  the  variables  k.  and  e.  The  default  value  used 
for  a£  in  the  AFVJAL/PORT  TEACH  version,  however,  was  1.21  and  all 
the  earlier  computations  had  been  performed  for  this  value.  We 
wanted  to  find  out  if  the  prediction  of  the  rear  stagnation  point 
is  influenced  by  a£.  For  a£  =  1.21,  the  normalized  forward  and 
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rear  stagnation  distances  (normalized  with  respect  to  the  center- 
body  diameter)  were  found  to  be  n.4565  and  1.2569.  Variation  of 
a£  from  0.90  to  1.40  in  steps  of  0.10  showed  that,  except  for 
j£  =  0.80  (when  both  fore  and  aft  stagnation  points  moved  farther 
downstream),  the  value  of  a£  did  not  have  much  effect  on  the  pre¬ 
dictions.  In  particular,  for  a£  =  1.3,  the  fore  and  aft  stagna¬ 
tion  points  were  at  the  normalized  distances  of  0.4508  and 
1.2575.  In  terms  of  the  ratio  of  the  stagnation  distances  and 
the  ratio  of  inlet  velocities,  this  result  corresponded  to  zp/z A 
=  0.3585  and  Wp/W ^  *  1.025  .  V7hile  the  results  for  a£  =  1.21  and 
a£  =  1.3  did  not  differ  significantly,  from  the  viewpoint  of  com¬ 
parisons  with  other  TEACH  predictions  all  the  subsequent  predic¬ 
tive  results  were  obtained  for  o£  =  1.3. 

3. 2. 2. 2  Influence  of  cy 

To  examine  the  influence  of  c u  a  set  of  computations  was 
carried  out  for  the  same  air  and  CO2  velocities  as  selected 
earlier.  The  value  of  a£  was  kept  at  1.21.  The  values  of 
c  y  tested  were  0.12,  0.15  ,  0.25  ,  0.35  ,  0.90  ,  and  1.20.  As  seen 
in  Figure  3,  cM  has  a  very  strong  influence  on  the  results. 
Increasing  the  value  of  c u  resulted  in  decreasing  both  the  fore 
and  aft  stagnation  distances.  Thus,  we  get  for  c y  =  0.12,  the 
values  of  0.3745,  1.1308,  and  0.3312,  and  for  cu  =  0.15,  the 
values  of  0.3308,  1.0509,  and  0.3148  as  the  normalized  forward 
stagnation  distance,  the  normalized  rear  stagnation  distance,  and 
their  ratios.  It  is  interesting  to  note  that  an  air  stagnation 
distance  of  roughly  one  centerbody  diameter  can  result  from  a 
cu  value  less  than  twice  the  accepted  value  of  0.09.  Since  an 
increase  in  c u  leads  to  an  increase  in  the  turbulent  eddy  viscos¬ 
ity,  it  would  seem  that  this  can  result  in  an  enhanced  turbulent 
mixing  rate  in  the  recirculation  region,  which  tends  to  mitigate 
the  overprediction  of  the  rear  stagnation  point.  No  doubt  this 
is  a  rather  crude  explanation,  in  view  of  the  fact  that  the 
influence  of  cu  is  pervasive  in  the  differential  equations  for 
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Figure  3.  Influence  of  cu  on  the  Centerline  Stagnation 
Points  (UCI  Combustor) . 
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axial  and  radial  momentum,  turbulent  kinetic  energy,  and  tur¬ 
bulent  dissipation. 

3. 2. 2. 3  Other  Effects 

All  the  above  computations  were  performed  with  the  41-X-34 
grid  and  with  the  exit-plane  boundary  of  the  computational  domain 
being  located  a  distance  of  6.2  cm  from  the  face  of  the  cen- 
terbody.  The  parameter  TURBIN  for  the  inlet  turbulent  kinetic 
energy  was  specified  as  0.03  and  the  parameter  ALAMDA  for  the 
inlet  turbulent  length  scale  was  specified  as  0.005.  A  limited 
number  of  sensitivity  tests  were  conducted  to  ascertain  the 
influence  of  grid  size  and  of  the  length-scale  parameter  ALAMDA. 
The  value  of  ALAMDA  employed  for  this  purpose  was  0.05  and  the 
coarse  axial  grid  intervals  corresponded  to  roughly  1.6  times  the 
previous  grid  intervals.  This  yielded  an  exit-plane  boundary 
location  of  10  cm.  These  tests  used  a£  =  1.3  and  cu  =  0.09.  The 
general  effect  of  the  ten-fold  increase  in  the  inlet  turbulence 
length  scale,  or  of  the  cruder  grid  size,  or  of  the  farther  exit 
boundary  was  to  move  both  the  forward  and  rear  stagnation  points 
slightly  upstream.  The  changes  were  relatively  insignificant, 
and  this  conclusion  conformed  to  the  results  of  the  sensitivity 
testing  reported  in  Paragraph  3.1.2  for  the  APL  configuration. 

3.3  SIMILARITY  CONSIDERATIONS 

The  preliminary  studies  of  the  centerbody  configuration 
indicated  that  the  numerical  predictions  resulted  in  qualita¬ 
tively  correct  flowfield  characteristics.  This  observation  was 
essentially  based  upon  the  axial  variation  of  the  centerline  mean 
axial  velocity  and  turbulence  kinetic  energy.  For  the  entire 
range  of  central  jet  flow  rates  from  the  annular  jet  dominant 
regime  (where  the  central  jet  is  completely  turned  back  toward 
the  centerbody)  to  the  central  jet  dominant  regime  (where  the 
central  jet  completely  penetrates  the  near-wake  recirculation 


region)  ,  the  numerical  results  were  consistent  with  the 
experimental  trends.  However,  as  emphasized  earlier,  the 
quantitative  results  of  the  centerline  stagnation  points  were  not 
satisfactory.  In  particular,  the  numerical  calculations 
underpredicted  the  forward  stagnation  points  and  overpredicted 
the  rear  stagnation  points.  Since  the  parameter  sensitivity 
tests  failed  to  explain  this  discrepancy  between  the  predictions 
and  measurements,  it  was  of  interest  to  ascertain  whether  the 
numerical  calculations  demonstrate  internal  consistency  with 
respect  to  certain  well-known  characteristics  of  free  turbulent 
shear  flows. 


Of  special  interest  was  the  radial  variation  of  the  mean 
axial  velocity  field.  Sufficient  experimental  evidence  exists 
for  the  similarity  of  mean  velocity  profiles  in  free  turbulent 
shear  flows  as  occurring  in  jets  and  wakes.  In  successive  axial 
stations  in  the  downstream  direction,  the  mean  axial  velocity 
profile  in  the  radial  direction  exhibits  the  more  or  less 
bell-shaped  curves  under  suitable  normalizations  of  the  velocity 
field  and  lateral  distance.  Since  the  similarity  of  mean  axial 
velocity  field  is  the  necessary  condition  for  correctly 
predicting  the  flowfield  behavior,  our  preliminary  numerical 
results  of  the  centerbody  flowfields  remained  to  be  checked  for 
this  aspect. 


3.3.1  Centerbody  Combustor  Configurations 


The  TEACH  code  predictions  of  the  radial  variation  of  the 
mean  axial  velocity  in  both  the  APL  and  UCI  configurations  were 
examined  from  the  viewpoint  of  obtaining  universal  profiles.  For 
this  purpose,  the  data  of  radial  coordinates,  r,  and  axial  velo¬ 
cities,  W(r),  were  normalized  following  the  suggestion  of 
Abramovich . ^2  The  normalized  radial  coordinate  is  defined  as  Ar 
=  ( r-r o . 5 ) / ( rg , g-rg . i )  and  the  normalized  axial  velocity  is 
defined  as  AW  =  LW(r)-Wmin  j/(wnax_wmin ' ’  Were,  Wmax  and  Wm^n  are 
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the  maximum  and  minimum  axial  velocities  at  any  given  axial  sta¬ 
tion.  The  values  of  rg  .  9,  rg.5  <  and  rg.^  correspond  to  the 
radial  locations  (at  any  given  axial  station)  where  (W-W^in)  is 
respectively  equal  to  90,  50,  and  10  percent  of  ('‘•max^'min )  (at 
the  same  axial  station). 

3. 3. 1.1  UCI  Configuration 

For  the  UCI  configuration,  experimental  conditions  in 
Reference  21  corresponded  to  a  reference  velocity  in  the  duct  of 
7.5  m/s  and  overall  equivalence  ratios  of  0.02  and  0.05.  For 
numerical  modeling,  two  different  annular  air  flow  rates  were 
considered,  corresponding  to  the  values  of  10.77  m/s  and  4.31  m/s 
for  the  annular  air  velocities.  Figures  4(a)  through  (g)  present 
the  results  for  the  air  velocity  of  10.77  m/s.  The  central-jet 
exit  velocity  for  COg  ranged  from  nearly  zero  to  30  m/s,  in  steps 
of  5  m/s.  For  all  cases  except  the  highest  COg  flow,  reverse  air 
flow  occurred,  giving  rise  to  two  stagnation  points  on  the  cen¬ 
terline.  As  anticipated  from  the  extrapolated  results  for  the 
centerline  peak  negative  axial  velocity  as  a  function  of  COg  exit 
velocity  which  gave  a  value  of  28  m/s  for  the  "breakthrough" 
velocity,  the  centerline  reverse  flow  was  completely  eliminated 
for  the  30  m/s  case. 

We  note  that  the  above  results  cover  the  entire  range  of 
conditions  from  the  annular  jet  dominant  case  to  the  central  jet 
dominant  case.  With  the  centerline  rear  stagnation  point  (when 
it  exists)  being  located  at  normalized  axial  distances  (z/D)  of 
1.2?  through  1.31,  the  radial  profiles  presented  correspond  to 
both  inside  and  outside  the  recirculation  region.  Finally,  in 
Figure  4(c)  through  (f)  (for  COg  exit  velocities  10  through 
25  m/s),  the  axial  stations  considered  are  upstream  and 
downstream  of  the  centerline  forward  stagnation  point.  An 
inspection  of  the  normalized  radial  profiles  of  the  mean  axial 
velocity  reveals  the  tendency  toward  a  universal  similarity 


profile.  The  solid  line  in  all  these  figures  is  the  result  for  a 
free  jet  (reproduced  from  Abramovich ) . 42 

It  must  be  noted  here  that  although  it  is  not  explicitly 
pointed  out  in  Abramovich,42  some  care  is  required  in  obtaining 
the  set  of  radii  rg.]_,  ^0.5/  and  rg.g.  It  appears  that  monoton¬ 
icity  is  a  necessary  condition  in  that  either  of  the  inequalities 

r0.1  >  r0 . 5  >  r0.9  or  r0.1  <  r0.5  <  r0.9  must  be  strictly  valid 
for  the  normalization  ( r-rg # 5 )/( rg ^ g-rg , 1 )  to  be  unambiguous.  It 
is  easy  to  see  that  the  former  inequality  applies  for  the  annular 
jet  near  the  confining  duct  and  the  latter  inequality  applies  for 
the  central  jet  near  the  centerline.  While  at  some  axial  sta¬ 
tions  and  under  certain  flow  conditions  a  single  set  of  radii  may 
suffice  all  the  way  from  the  centerline  to  the  duct  wall,  in 
general  it  has  been  found  that  upstream  of  the  rear  stagnation 
point  both  inequalities  occur.  When  two  such  sets  of  radii 
exist,  the  resulting  normalizations  conform  properly  to  the  uni¬ 
versal  profile  in  their  respective  regions  and  deviate  from  it 
in  the  other.  For  example,  one  set  of  radii  is  appropriate  for 
-1.25  <  Ar  <  0  and  the  other  for  0  <  Ar  <1.25.  The  results  pre¬ 
sented  here  are  the  composites  of  two  such  normalizations  with 
the  deviant  portions  of  the  curve  being  deleted  in  the  figures. 
Although  the  implications  of  these  observations  are  not  clear,  it 
seems  that  jet-like  and  wake-like  behaviors  occur  in  the 
appropriate  regions. 

Furthermore,  for  the  free  jet  issuing  into  a  quiescent 
region,  Wmax  occurs  at  the  centerline  and  Wm^n  (=  0)  occurs 
asymptotically  at  the  "edge"  of  the  jet.  Thus,  both  AW  and  Ar 
are  unambiguous.  The  juxtaposition,  however,  of  the  confining 
duct,  annular  jet,  centerbody-wake  region,  and  central  get  does 
introduce  an  essential  element  of  nonuniqueness.  For  example, 
when  the  central  jet  still  retains  its  forward  momentum,  the 
centerline  axial  velocity  is  positive  and  is  also  greater  than 
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the  axial  velocities  at  neighboring  off-centerline  radial 
locations.  Thus,  the  centerline  velocities  represent  Wmax 
those  axial  stations.  Downstream  of  the  forward  stagnation 
point,  the  centerline  axial  velocities  are  negative  and 
correspond  to  Wmj_n  for  the  axial  locations  between  the  forward 
and  rear  stagnation  points.  It  is  also  easy  to  see  that  upstream 
of  the  forward  stagnation  point  there  are  of f-centerl ine  radial 
locations  where  the  axial  velocities  are  negative  and  represent 
W^in  at  those  axial  locations.  Similar  off-centerline  obser¬ 
vations  can  be  made  when  the  central  jet  has  eliminated  the 
reverse  flow  on  the  centerline.  These  differences  account  for 
the  discrepancies  between  the  normalized  curves  for  the  free  jet 
and  the  centerbody  flowfield. 

The  above  comments  are  further  confirmed  when  we  examine  the 
behavior  of  the  universal  curve  as  a  function  of  the  exit 
velocity  of  the  CO2  jet.  With  increasing  central  jet  velocity, 
the  tendency  of  the  predicted  data  to  show  greater  conformity 
with  the  free  jet  curve  is  unmistakable.  Also,  we  notice  the 
velocity  "overshoot"  for  0.25  <  Ar  <  0.50  for  the  different  axial 
locations  (the  closer  the  axial .  location  to  the  centerbody,  the 
larger  the  discrepancy).  An  examination  of  the  raw  data  shows 
that  this  region  denotes  the  entrainment  of  the  CO2  jet  by  the 
shear  layer  of  the  annular  jet.  Clearly  this  effect  should 
diminish  as  the  strength  of  the  central  jet  increases.  This  is 
evident  from  Figure  4(f)  and  (g). 

In  Figures  5(a)  through  (f),  the  predicted  results  for  the 
air  velocity  of  4.31  m/s  are  shown  for  the  central-jet  exit 
velocity  for  CO2  ranging  from  nearly  zero  to  15  m/s.  Since  the 
critical  exit  velocity  for  "breaking  through"  the  recirculation 
reaion  is  found  to  be  12  m/s  by  extrapolation,  the  flowfield  for 
15  m/s  CO2  exit  velocity  represents  the  central  jet  dominant 


•51 


-0.6  -0.4 


.6  0. 


Figure 


-0.8 


-0.2  0  0.2  0.4  0 

Ar 


(d) 


5.  Mean  Axial  Velocity  Preditions  for  the  (UCI) 
Centerbody  Configuration. 


regime  wherein  the  reverse  flow  along  the  centerline  is 
eliminated.  The  normalized  profiles  exhibit  similarity  as  before 
and  these  data  also  confirm  that  the  universality  of  the  velocity 
profile  holds  at  different  annular  air  velocities  as  well. 

3. 3. 1.2  APL  Configuration 

The  predicted  results  for  the  APL  configuration  are  shown  in 
Figures  6(a)  through  (c).  Here,  the  flowfields  considered  are 
those  for  the  annular  air  mass  flow  of  2  kg/s  and  central  CO2 
mass  flows  of  4,  8,  and  16  kg/hr.  Although  TEACH  computations 
were  completed  for  other  annular  air  mass  flows  of  0.07,  0.5,  1, 
and  3  kg/s  and  for  a  large  range  of  CO2  mass  flows  in  each  case, 
Abramovich -type  normalizations  were  not  done  for  all  these  cases. 
However,  from  the  smaller  subset  seen  in  Figure  6,  it  is  clear 
that  combustor  scaling  preserves  the  similarity  of  the  radial 
profile  of  the  mean  axial  velocities. 

The  foregoing  observations  based  upon  TEACH  code  predictions 
are  strengthened  by  an  examination  of  experimental  data  of  the 
profile  measurements  in  the  APL  configuration.  The  normalization 
of  all  the  available  data-'-2,4^  for  annular  air  mass  flow  of 
2  kg/s  is  seen  in  Figure  7.  It  is  clear  from  these  data  that  the 
predicted  behavior  is  real  in  the  centerbody  configuration. 

3.3.2  Other  Flowfield  Configurations 

In  view  of  the  observed  universal  nature  of  the  radial 
variation  of  the  axial  velocity  in  the  centerbody  configuration, 
it  is  of  interest  to  examine  other  flowfields  investigated  in  the 
literature.  Figure  8,  which  shows  the  configuration  and  the 
normalized  results  of  Abramovich  and  Vafin,42  forms  the  basis  of 
our  line  of  inquiry.  This  configuration  is  quite  similar  to  that 
of  the  centerbody  combustor.  It  differs  only  in  the  absence  of 
the  central  jet  and  in  the  much  shorter  length-to-diameter  ratio 
of  the  centerbody.  The  previously  noted  departure  of 
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experimental  data  from  the  free  jet  behavior  is  clearly  observed 
here.  Furthermore,  one  notices  the  large  acceleration  of  the 
annular  flow  and  the  consequent  reduction  of  the  boundary-layer 
growth  in  this  configuration.  It  appears  that  from  the  results 
of  our  centerbody  configuration,  which  has  a  greater  propensity 
for  boundary-layer  growth,  the  similarity  of  the  mean  velocity 
profiles  does  not  depend  strongly  on  the  presence  of  the 
centerbody  boundary  layer. 

The  flowfields  discussed  thus  far  represent  confined  flows 
because  of  the  presence  of  the  outer  duct  wall.  An  example  of  an 
unconfined  flowfield  is  that  downstream  of  an  annular  nozzle, 
investigated  by  Ko  and  Chan. While  the  complete  details  of 
their  configuration  are  not  available,  the  annular  jet  is  formed 
by  the  flow  of  a  uniform  stream  over  one  end  of  a  long  cylinder, 
with  its  axis  aligned  parallel  to  the  direction  of  the  stream  and 
located  concentr ically  inside  a  converging  duct.  The  trailing 
end  of  the  cylinder  is  made  flush  with  the  exit  plane  of  the  duct 
(see  the  upper  part  of  Figure  9).  The  outer  diameter,  0o,  of  the 
annular  jet  is  6.2  cm  and  the  inner  diameter,  Dj_ ,  (which  is  the 
diameter  of  the  cylinder)  is  2.8  cm,  giving  a  blockage  ratio 
(=  1  -  D?/D^)  of  0.80  (the  value  quoted  in  the  paper  is  0.78). 

For  the  APL  conf iguration ,  the  blockage  ratio  is  0.70. 

The  hot-wire  measurements  of  the  axial  velocity  as  a 
function  of  the  radial  position  at  different  axial  locations  are 
shown  in  the  normalized  coordinates  in  the  lower  part  of  Figure 
9.  Ko  and  Chan  superimposed  the  confined  annular  jet  results  of 
Abramovich  and  Vafin^  as  well  as  the  results  of  the  single  jet 
and  demonstrated  the  existence  of  similarity  of  the  mean  velocity 
profiles  for  their  unconfined  configuration.  We  note  that  the 
farthest  axial  station  of  their  measurements  is  at  1.1  cylinder 
diameter,  just  short  of  the  centerline  reattachment  point 
(located  at  1.11  cylinder  diameter).  Our  predictions  for  the 
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centerbody  configurations,  however,  extend  to  more  than  two 
diameters  of  the  centerbody  (thus,  farther  downstream  of  tne 
reattachment  point)  and  still  retain  the  similarity  of  the  axial 
velocity  field. 

Although  the  upstream  details  of  the  forebody  and  the  outer 
converging  nozzle  are  not  provided  by  Ko  and  Chan,  we  may 
conclude  that  their  annular  jet  stream  is  nearly  axial  at  the 
exit  plane.  While  the  boundary-layer  growth  on  the  inner  wall  of 
the  nozzle  is  expected  to  be  minimal,  there  can  be  appreciable 
boundary-layer  development  on  the  cylinder,  depending  on  how  long 
it  is.  An  exit  configuration  exhibiting  little  boundary-layer 
effects  but  which  introduces  a  radial  velocity  component  at  the 
exit  plane  is  the  annular  jet  of  Durao  and  VIhitelaw,45 
illustrated  in  the  upper  part  of  Figure  10(a).  This  is  again  an 
unconfined  flowfield,  downstream  of  a  disk  concentrically  located 
inside  a  converging  nozzle.  Laser  Doppler  velocimetry  was 
employed  to  furnish  the  axial  and  radial  velocity  components 
downstream  of  the  disk.  The  exit-plane  measurements  revealed  a 
significant  radial  velocity  component.  We  have  presented  their 
raw  data  of  the  radial  variation  of  the  mean  axial  velocity  in 
the  normalized  form,  as  seen  in  the  lower  part  of  figure  10(a). 
The  results  shown  here  pertain  to  the  disk  of  14.2  mm  diameter 
(giving  a  blockage  ratio  of  0.50)  and  an  exit  axial  velocity  of 
26.3  m/s.  For  these  conditions,  the  centerline  reattachment 
point  is  located  at  one  disk  diameter  downstream.  The  normalized 
results  show  that  the  velocity  profiles  in  the  near  field 
(upstream  of  the  stagnation  point)  do  not  exhibit  similarity. 

The  deletion  of  the  data  corresponding  to  the  first  two  axial 
stations  (i.e.,  at  0.14  and  0.42  disk  diameters)  as  illustrated 
in  Figure  10(b)  clarifies  this  point.  We  believe  that  the 
presence  of  significant  outward  radial  velocity  component  in  the 
annular  stream  very  close  to  the  exit  plane  distinguishes  this 
configuration  from  all  the  other  configurations  discussed  earlier 


72 


Annular-jet  configuration.  Dimensions  in  mm.  —  39-5,  26-8,  16-9  or  9-4  m/s; 
D  =  20-0  mm;  d  =  14-2,  12-3  or  8-9  mm;  disks  manufactured  with  sharp  edges. 
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Figure  10  (b) .  Mean  Axial  Velocity  Measurements  for  the 
Unconfined  Annular  Jet  Around  a  Disk. 


and  thereby  contributes  to  the  marked  departure  from  similarity. 
Farther  downstream,  where  the  radial  velocity  components  are 
directed  inward  (toward  the  centerline),  the  flowfield  begins  to 
resemble  all  the  near-wake  flows  considered  earlier  (presumably, 
the  downstream  flow  development  no  longer  retains  the  memory  of 
the  initial  profile  at  the  exit  plane). 

The  foregoing  lends  support  to  our  viewpoint  that  all  the 
flowfields  investigated  so  far  belong  to  one  class  of  turbulent 
flows  which  obeys  certain  similarity  considerations  for  the  mean 
axial  velocity.  The  composite  of  all  the  data  (predictions  and 
measurements)  for  all  the  configurations  displayed  in  Figure  11 
vividly  demonstrates  this  point. 

3.3.3  Species  Concentration  Fields 

Previous  fluid  dynamic  studies  have  shown  that  when 
temperature  and  species  concentration  effects  are  present, 
suitably  normalized  temperature  and  concentration  variables  also 
exhibit  similarity  (see  Ref.  42,  for  example).  In  our  isothermal 
predictive  modeling  of  the  turbulent  mixing  of  air  and  CC>2/  the 
results  of  the  radial  variation  of  CO2  mass  fractions  at 
different  axial  locations  are  available.  Abramovich-type 
normalizations  for  the  mass  fraction  profiles,  however,  have 
failed  to  yield  the  universal  profiles.  Profiles  of  variables 
such  as  the  mole  fractions  of  CO2  or  the  concentration  (in 
moles/unit  volume)  of  CO2  have  also  been  unable  to  conform  to 
similarity  considerations. 

We  found  these  results  somewhat  puzzling.  Having 
established  the  similarity  of  axial  velocity  fields,  it  was 
difficult  to  accept  that  the  species  field  behaves  differently. 

A  possible  source  of  this  discrepancy  is  the  use  of  unity  Schmidt 
number  for  determining  the  turbulence  exchange  coefficient  for 
the  species  field  from  the  values  of  turbulent  eddy  viscosity  in 
the  predictive  calculations .  However,  this  effect  is  expected  to 
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Figure  11.  A  Composite  of  all  the  Predicted  and  Measured 
Results. 
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be  generally  insignificant.  A  more  serious  pitfall  may  be  the 
rather  small  amount  of  CO2  present  (compare  a  2  kg/s  air  mass 
flow  with  an  8  kg/hr  CO2  mass  flow).  The  numerical  results 
indicated  that  the  conservation  of  CO2  mass  was  not  satisfied 
over  the  duct  cross  section  at  each  axial  location — the  CO2  mass 
flow  was  underestimated  by  as  much  as  50  percent,  especially  near 
the  forward  stagnation  point.  While  the  velocity  field  remained 
largely  unaffected  (as  evidenced  by  the  similarity)  by  the 
failure  of  the  CO2  mass  conservation  (in  view  of  the  trace  amount 
of  CO2  overall),  this  was  not  likely  to  be  true  of  the  radial 
variation  of  CO2  concentrations.  Clearly  this  aspect  needed 
further  study.  Whether  the  lack  of  similarity  for  CO2 
concentrations  was  due  to  the  inadequacy  of  numerical  modeling 
may  perhaps  be  ascertained  by  examining  the  APL  experimental  data 
(of  probe  sampling  of  CO2)  for  Abramovich-type  normalizations. 
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SECTION  4 

REFINEMENTS  AND  SELECTED  RESULTS 


Although  the  computational  results  reported  in  the  previous 
section  demonstrated  that  the  TEACH-type  numerics  and  the  k-e 
turbulence  model  possessed  the  capability  to  predict  physically 
correct  flowfields  in  the  centerbody  combustor  configuration,  the 
comparison  with  the  experimental  data  found  these  predictions 
wanting.  Further  computational  investigations  were  necessary  for 
obtaining  improved  predictions.  Several  areas  of  refinement  were 
examined  and  these  are  discussed  in  this  section. 

4.1  INFLUENCE  OF  DIFFERENCING  SCHEMES 

Some  TEACH  code  calculations  were  carried  out  to  study  the 
effects  of  different  differencing  schemes  on  the  centerline 
locations  of  stagnation  points.  It  had  been  noted  earlier  that 
the  TEACH  code  results  underpredicted  the  forward  stagnation 
point  and  overpredicted  the  rear  stagnation  point.  While  the 
latter  is  expected  to  be  affected  by  the  strong  streamline 
curvature  in  the  vicinity  of  the  rear  stagnation  point  and  the 
failure  of  the  standard  k-e  model  to  include  this  effect,  the 
amount  of  numerical  diffusion  inherent  in  various  differencing 
schemes  appears  to  affect  significantly  the  former. 

In  the  TEACH  code,  the  convective  terms  have  been  discre¬ 
tized  through  the  "hybrid”  upwind  differencing  scheme  (see 
Paragraph  2.4.2).  We  completed  some  calculations  by  replacing 
this  scheme  with  the  power-law  differencing  scheme  (see  Paragraph 
2.4.3).  For  the  case  of  2  kg/s  air  flow  and  3  kg/hr  CO2  flow, 
the  centerline  locations  of  the  forward  and  rear  stagnation 
points  (in  terms  of  the  centerbody  diameter)  were  as  follows. 
Forward:  0.-101  for  "hyorid"  upwind  with  fine  grid,  and  0.423  and 

0.44  for  power  law  with  fine  and  coarse  grids  respectively. 

Rear:  1.347  for  "hybrid"  upwind  with  fine  grid,  and  1.340  and 
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1.314  for  power-law  with  fine  and  coarse  grids  respectively.  We 
note  that  the  power-law  scheme  results  in  7  to  10  percent 
increase  in  the  forward  stagnation  point  distance  and  0.5  to  2.5 
percent  decrease  in  the  rear  stagnation  point  distance.  3otn 
these  changes  are  in  the  desired  direction  but,  as  pointed  out 
earlier,  the  rear  stagnation  point  is  less  sensitive  to  the 
change  in  the  differencing  scheme. 

4.2  C02  MASS  CONSERVATION 

Additional  computations  were  carried  out  to  examine  the  mass 
entrainment  of  one  jet  by  the  other.  For  the  ducted  flows 
exemplified  by  the  centerbody  combustor  configuration,  the  con¬ 
cept  of  entrainment  is  different  from  that  of  unconfined  (free) 
jet  flows  where  additional  mass  is  entrained  from  the  ambient 
fluid.  Here,  the  confining  boundary  implies  that  there  is  no 
additional  mass  being  entrained,  but  due  to  turbulent  convection 
and  diffusion,  there  is  a  redistribution  of  the  air  and  C02  mass 
fluxes  at  different  axial  and  radial  locations.  This  is  readily 
seen  in  Figure  12.  We  note  that  when  the  annular  air  flow 
dominates  the  flowfield  (encountered  typically  at  air  flow  rates 
of  2  kg/s  and  C02  flow  rates  of  4  or  8  kg/hr),  the  central  C02 
flow  is  brought  to  rest  on  the  centerline  by  the  reverse  flow  of 
air,  carried  radially  outward,  and  entrained  by  the  inner  region 
of  the  annular  air. 

Before  these  local  redistributions  can  be  quantified,  it  was 
necessary  to  assure  ourselves  that  the  global  mass  conservation 
was  satisfied.  Unfortunately,  our  integrated  mass-flow  results 
revealed  that  while  the  annular  air  mass  flow  is  conserved  at  all 
axial  locations,  the  central  C02  mass  flow  failed  to  satisfy  mass 
conservation.  At  axial  locations  in  the  vicinity  of  the  cen¬ 
terline  forward  stagnation  point,  the  integrated  C02  mass  flow 
was  underestimated  by  as  much  as  50  percent.  Our  preliminary 
results  indicated  that  this  discrepancy  arose  in  both  "hybrid" - 
upwind  and  power-law  differencing  schemes. 
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An  examination  of  the  code  showed  that  the  mass  conservation 
was  satisfied  only  for  the  air./CC>2  mixture  through  the  overall 
continuity  equation  in  the  so-called  SIMPLE  procedure^  for 
pressure  correction.  Note  that  for  CO2  mass  flows  of  3  kg/hr  and 
less,  the  ratio  of  CO2  mass  flow  to  air  mass  flow  is  of  the  order 
of  10 “3  when  the  air  flow  is  2  kg/s.  Thus,  it  appeared  that  the 
failure  to  satisfy  COg  mass  conservation  might  be  traceable  to 
the  insignificant  contribution  made  by  the  CO2  jet  to  the  total 
mass  flow.  Figure  13  shows  the  integrated  axial  mass  flow  of  air 
and  CO2  (per  unit  radian  of  the  cross  section)  at  different  axial 
locations.  Our  first  set  of  TEACH  code  computations  employed  the 
value  of  0.05  for  the  inlet  turbulence-length-scale  parameter  X 
and  the  computational  grid  denoted  A.  The  nodal  distribution  of 
this  grid  is  shown  in  Figure  1.  Figure  13  shows  that  air  mass 
flow  satisfies  the  conservation  requirement  very  well.  The  COg 
mass  flow,  however,  increasingly  diverges  from  the  conserved 
value  and  is  off  by  nearly  20  percent  (of  the  inlet  value)  at  the 
exit  boundary  of  the  computational  domain.  At  the  cross  section 
corresponding  to  the  centerline  forward  stagnation  point,  the 
integrated  CO2  mass  flow  attains  a  minimum. 

It  seemed  to  '  -hat  the  observed  discrepancy  in  COg  mass 
conservation  may  be  due,  at  least  in  part,  to  poor  spatial 
resolution  of  the  computational  grid.  In  Figure  1  we  see  that 
the  grid  A  is  characterized  by  densely  populated  grid  nodes  near 
the  centerline  and  centerbody  face  and  rather  sparsely  populated 
nodes  near  the  duct  boundary  and  exit.  In  particular,  the 
annular  region  has  only  eight  interior  grid  points  in  the  radial 
direction.  Although  the  spatial  resolution  is  quite  good  toward 
the  bottom  left  of  the  computational  domain,  which  is  normally 
the  region  of  interest  for  the  CO2  jet,  under  the  annular 
dominant  regime  considered  here  the  CO2  jet  is  turned  into  the 
annular  region.  Thus,  the  calculations  for  the  CO2  mass  fraction 
and  the  integrated  mass  flow  are  greatly  affected  by  the  poor 
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spatial  resolution  in  this  region,  especially  in  view  of  the 
rather  small  value  of  CO2  mass  flux.  Therefore,  a  distribution 
of  grid  nodes  showing  better  spatial  resolution  in  the  annular 
region  can  be  expected  to  result  in  less  discrepancy  with  the 
conservation  requirement.  This  conjecture  also  seemed  plausible 
in  view  of  the  calculations  for  the  central  jet  dominant  regime 
(wherein  the  CO2  jet  "breaks  through"  the  bluff-body 
recirculation  region)  which  satisfied  the  mass  conservation 
requirement . 

To  test  this  hypothesis,  we  completed  calculations  with  the 
grid  denoted  B  shown  in  Figure  14.  By  a  redistribution  of  the 
grid  nodes  in  the  radial  direction,  grid  3  resulted  in  thirteen 
nodes  in  the  interior  of  the  annular  region.  The  results  for 
this  case  (for  A  =  0.05)  seen  in  Figure  13  indicate  that  our 
conjecture  is  essentially  correct.  The  conservation  of  CO2  mass 
flow  is  much  better  everywhere,  except  near  the  forward  stagna¬ 
tion  point.  That  the  spatial  resolution  of  the  computational 
grid  may  not  be  the  only  factor  contributing  to  the  question  of 
mass  conservation  was  confirmed  by  the  results  shown  in  Figure  13 
for  A  =  0.3333  (the  rationale  for  specifying  this  value  is  pre¬ 
sented  subsequently).  However,  unlike  the  grid  distribution 
whose  influence  is  largely  numerical,  the  increase  in  A  implies 
large  inlet  turbulence  length  scale  which  results  in  enhanced 
turbulence  activity  in  the  flowfield.  This  is  expected  to 
augment  the  turbulent  diffusion  near  the  forward  stagnation 
point,  thereby  reducing  the  discrepancy  in  mass  conservation  in 
that  vicinity. 

From  the  above  discussion  it  would  seem  that  the  predictive 
results  would  benefit  by  increased  spatial  resolution  and  thus  by 
the  use  of  a  greater  number  of  grid  nodes  than  have  been  employed 
in  our  numerical  experiments.  The  observed  minimum  at  the  for¬ 
ward  stagnation  point  and  the  peak  between  the  jet  exit  and  the 
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Figure  14.  An  Arbitrary  Modification  of  the  Grid  with  Increased 
Spatial  Resolution  in  the  Annular  Region  (GRID  B) . 
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forward  stagnation  point,  however,  do  not  seen  to  be  entirely 
related  to  the  question  of  spatial  resolution. 

To  explore  this  aspect  further,  we  have  shown  in  Figure  15 
the  relative  contribution  of  the  diffusive  and  convective  fluxes 
to  the  integrated  CO2  mass  flow  at  different  axial  locations. 

This  shows  that  at  the  cross  section  near  the  forward  stagnation 
point,  the  influence  of  diffusion  is  the  highest,  and  that  while 
the  CC>2  mass  flow  in  the  axial  direction  is  largely  due  to  con¬ 
vective  transport  in  a  large  part  of  the  flowfield,  there  exist 
regions  where  diffusive  contribution  is  no  longer  negligible. 
Neglect  of  diffusive  flux  can  result  in  a  large  discrepancy  in 
the  COj  mass  conservation,  and  was  seen  to  result  in  as  much  as 
50  percent  underprediction  of  the  integrated  CO2  mass  flow  near 
the  forward  stagnation  point,  as  mentioned  previously.  The 
results  in  Figure  13,  however,  include  the  diffusive  contribution 
to  the  total  CO2  mass  flow.  Thus,  the  peculiar  trend  observed 
upstream  of  the  forward  stagnation  point  remains  unclear,  except 
for  possible  spurious  effects  arising  from  numerical  diffusion. 

4.3  INFLUENCE  OF  INLET  TURBULENCE  LENGTH  SCALES  REVISITED 

A  number  of  recent  turbulent  recirculating  flowfield 
calculations  have  underpredicted  the  extent  of  the  recirculation 
region,  the  degree  of  underprediction  depending  on  the 
configuration,  and  the  turbulence  model  employed.  Our  isothermal 
flowfield  modeling  of  the  cer.terbody  configuration  with  the 
version  of  the  TEACH  code  available  to  us,  however,  has  resulted 
in  an  overprediction  of  the  centerline  location  of  the  reattach¬ 
ment  point  due  to  the  annular  jet  around  the  centerbody.  At  the 
same  time,  a  second  stagnation  point  located  farther  upstream, 
due  to  the  introduction  of  a  weak  central  jet,  has  been  generally 
underpredicted.  The  numerical  results  reported  by  Sturgess  and 
Syed^  also  exhibit  this  overprediction  of  the  rear  stagnation 
point  and  underprediction  of  the  forward  stagnation  point.  3ut 
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Figure  15.  The  Relative  Importance  of  Diffusion  to  the  Total 
CO-,  Mass  Flow  Rate  at  Different  Axial  Locations. 
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the  degree  of  overprediction  of  the  rear  stagnation  point 
location  was  much  less  than  that  seen  in  our  results.  Since  our 
preliminary  sensitivity  tests  (discussed  in  Paragraphs  3.1.2  and 
3.2.2)  with  respect  to  several  parameters  of  the  numerical  model 
had  not  explained  the  source  of  the  discrepancy,  it  was  of 
interest  to  examine  this  question  further. 

Our  subsequent  examination  of  the  code  and  the  k- e  model 
therein  revealed  that  the  answer  to  the  above  difficulty  lay  in 
the  specification  of  the  inlet  turbulence  length  scales.  This 
subject  has  been  discussed  at  length  in  Paragraphs  2.2  and  2.4 
and  it  was  essential  that  the  distinction  between  the  length 
scales  anf?  l2  lsee  ^qs*  (7b)  and  (9b)  respectively  be  taken 
into  account  in  the  proper  specification  of  the  inlet  length 
scale . 

It  was  reported  by  Sturgess  and  Syed^  that  the  rear  stagna¬ 
tion  point  location  on  the  centerline  is  very  sensitive  to  the 
inlet  length  scale  and  hence  the  value  of  X.  The  value  of  X  used 
in  their  calculations  was  0.03.  Our  earlier  sensitivity  tests 
examined  the  range  of  0.005  to  0.05  (see  Tables  7  and  3)  for  the 
value  of  X  and  discerned  very  little  influence  on  the  stagnation 
point  location.  However,  as  noted  previously,  different  for¬ 
mulations  (  or  1 2)  have  been  used  in  these  studies.  Since  £]_  = 

%2/<z  u  [and  2-2  =  0*03  (R2_Ri)  ],  our  previous  parametric  range  of 
0.005  through  0.05  was  very  much  smaller  than  a  value  of  X  = 
0.3333  which  arises  from  the  requirement  for  consistency  ti.e., 
i]_  =  X  (R2-Ri>  =  2.2/0'09  *  °*°3  (R2_r1>/0*  09  =  0*  3333  (R2-Ri)  ]. 
Accordingly,  the  problem  of  centerline  axial  velocity  charac¬ 
teristics  and  stagnation-point  locations  was  reexamined,  in 
anticipation  that  the  TEACH  predictions  would  show  better 
comparison  with  the  experimental  data.  30 

Numerical  experiments  were  conducted  with  different  values 
of  X,  viz.,  0.05,  0.3333,  and  0.5556.  Note  that  these 


corresponded  respectively  to  the  arbitrarily  small  value  used  in 
our  earlier  studies,  the  "standard"  value  of  Sturgess  and  Syed® 
and  the  value  used  by  Leschziner  and  Rod i. 33  Furthermore,  with 
c u  =  0.09,  X  =  0.3333  corresponds  to  the  length  scale  of  three 
percent  of  the  reference  length  and  X  =  0.5556  corresponds  to 
five  percent  of  it  (0.05/0.09  =  0.5556). 

Figure  16  shows  the  results  of  the  centerline  forward  and 

rear  stagnation  points  for  a  fixed  annular  air  flow  of  2  kg/s  and 

a  number  of  central  CO2  flow  rates.  Identical  values  of  ALAMDA 

and  TURBIN  (0.03)  in  our  numerical  experiments  and  those  of 

/* 

Sturgess  and  Syed0  facilitate  the  comparison.  The  grid  systems  A 
and  B  in  our  experiments  are  those  denoted  in  Figures  1  and  14  (B 
displaying  a  more  dense  grid  in  the  annular  region  than  A).  Vie 
note  that  our  results  with  the  grid  system  A  show  good  agreement 
with  those  of  Sturgess  and  Syed.®  It  would  seem  that  their  grid 
system  is  very  similar  to  ours.  However,  it  should  be  noted  that 
their  results  are  based  upon  the  computational  domain  with  the 
exit-plane  boundary  located  at  0.8  m  from  the  centerbody  face.® 
However,  this  boundary  is  only  0.3  m  away  in  our  calculations . 
Therefore,  it  appears  that  the  computed  results  are  not 
significantly  affected  over  this  range  of  exit-plane  locations, 
so  long  as  the  exit  boundary  is  located  sufficiently  downstream 
of  the  rear  stagnation  point.  To  the  extent  that  numerical 
predictions  show  excellent  agreement  with  each  other,  the  degrees 
of  their  underprediction  of  the  forward  stagnation  points  and  the 
overpred iction  of  the  rear  stagnation  points  with  respect  to  the 
measured  data^-O  are  similar. 

Our  results  with  the  two  different  grid  systems  indicate  the 
appreciable  grid  sensitivity,  especially  for  the  rear  stagnation 
points.  Although  the  choice  of  grid  3  resulted  in  improved  CO2 
mass  conservation  (see  Figure  13),  its  predictions  of  the  rear 
stagnation  points  are  worse  than  those  of  grid  A.  In  view  of  the 
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fixed  number  of  grid  nodes  employed  in  these  two  grid  systems,  a 
mere  redistribution  of  the  grid  points  over  different  regions  of 
the  flowfield  can  only  result  in  different  aspects  of  the 
flowfield  predictions  displaying  varying  degrees  of  "success. " 
What  is  clearly  established  here  is  that  the  numerical  modeling 
could  do  better  with  increased  spatial  resolution  everywhere  in 
the  computational  domain.  That  the  predictions  of  the  forward 
stagnation  points  by  the  two  grid  systems  do  not  significantly 
differ  from  each  other  is  merely  a  consequence  of  the  fact  that 
the  spatial  resolution  in  the  vicinity  of  the  central  jet  is 
essentially  the  same  for  both  grid  systems. 

Figure  17  presents  the  centerline  stagnation  point  results 
for  different  values  of  a.  The  grid  system  B  is  chosen  in  this 
comparison  in  anticipation  of  the  subsequent  discussion  on  the 
curvature  correction  (whose  performance  is  better  evaluated  with 
the  poorer  predictions  of  grid  B)  .  We  observe  that  as  X 
increases,  the  rear  stagnation  point  moves  upstream  signifi¬ 
cantly.  This  behavior  has  been  noticed  by  Sturgess  and  Syed^ 
also.  The  effect  on  the  forward  stagnation  point  is  not  that 
strong,  especially  since  the  curvature  correction  with  A  =  0.5556 
tends  to  result  in  an  opposite  effect  (although  in  better 
agreement  with  the  experimental  data^-®).  Of  greater  significance 
is  the  dramatic  reduction  in  the  overprediction  of  the  rear 
stagnation  point  as  A  increases  from  0.05  (a  value  used  in  our 
earlier  studies)  to  0.3333.  Thus,  a  major  discrepancy  noted 
earlier  between  our  TEACH  predictions  and  those  reported 
elsewhere^  has  been  resolved  to  our  satisfaction. 

We  note  here  that  the  influence  of  the  inlet  length  scale  on 
the  forward  and  rear  stagnation  points  recalls  our  earlier 
numerical  experiments  (discussed  in  Paragraph  3. 2. 2. 2)  to 
determine  the  effect  of  cu.  These  calculations  were  carried  out 
for  the  UCI  configuration  which  is  roughly  a  1/5-scale  model  of 
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Figure  17.  Influence  of  Inlet  Turbulence  Length  Scale  and  Curvature  Correction 
on  the  Centerline  Stagnation  Points. 


the  APL  configuration.  These  results  are  shown  in  Figure  3.  We 
note  that  although  the  values  of  X  and  aE  employed  in  this 
experiment  are  different  from  those  in  Figure  17,  the  trend  is 
similar  to  the  length-scale  experiment,  in  that  both  stagnation 
points  move  upstream,  with  the  rear  stagnation  point  more 
strongly  affected  than  the  forward  stagnation  point.  However,  it 
must  be  pointed  out  that  in  the  k- e  model,  the  equilibrium- 
turbulence  assumption  implies  that  the  two  constants  cu  and  cD 
vary  together,  so  that  cu  •  Cp  =  0.09  (see  Reference  28).  Our 
computational  experiment  varied  only  cy. 

It  seems  that  larger  values  of  c u  or  larger  values  of 
inlet-length  scale  imply  higher  values  of  turbulence  eddy 
viscosity.  Therefore,  the  results  obtained  with  the  variation  of 
c u  or  X  (Figures  3  and  17)  suggest  that  a  higher  level  of 
turbulence  activity  results  in  shorter  recirculation  length. 

There  seems  to  be  some  experimental  evidence  to  support  this 
conclusion.  For  example,  an  experimental  comparison  46  of 
cold-flow  and  combusting-f low  recirculation  lengths  revealed  that 
the  turbulence  level  was  much  lower  and  the  recirculation  length 
was  much  larger  in  combusting  flows  than  in  cold  flows.  The 
centerline  measurements-'-0  of  the  axial  components  of  the  mean 
velocity  and  turbulent  intensity  in  isothermal  and  combusting 
flows  on  the  centerbody  combustor  configuration,  however,  do  not 
demonstrate  this  behavior.  While  the  recirculation  lengths  are 
indeed  larger  in  combusting  flows  than  in  cold  flows,  there  is  no 
clear  evidence  that  the  combusting  flows  indicate  decreased 
turbulence  activity  in  comparison  with  the  cold  flow. 

4.4  STREAMLINE  CURVATURE  EFFECTS  IN  TURBULENCE  MODELING 

That  the  streamline  curvature  has  a  strong  influence  on  the 
shear-flow  turbulence  is  well  known  in  the  literature.  Indeed, 
this  was  the  topic  in  the  comprehensive  monograph  by  Bradshaw. 32 
Since  the  size  of  the  recirculation  region  appeared  to  depend 
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strongly  on  the  turbulence  activity  (as  indicated  in  the  previous 
section)  in  the  curved  shear  layers  bordering  this  region,  and 
since  the  "standard"  k- s  model  (discussed  in  Paragraph  2.2)  does 
not  account  for  the  curvature  effects,  it  was  considered 
worthwhile  to  reexamine  the  discrepancies  between  predicted  and 
measured  results  in  the  light  of  streamline-curvature  corrections 
implemented  in  the  code  (as  outlined  in  Paragraph  2.3). 

Numerical  calculations  were  carried  out  by  introducing  a 
curvature-dependent  cu  (see  Paragraph  2.3  for  details)  into  the 
"standard"  k- e  model,  along  the  lines  suggested  by  Leschziner  i 
Rodi.-^  we  emphasized  earlier  the  ad  hoc  nature  of  these  corr 
tions.  Furthermore,  our  numerical  experiments  are  not  complet* 
and  only  a  limited  parametric  range  of  centerbodv-f lowf ield  co 
ditions  has  been  considered  so  far.  However,  the  results 
reported  here  do  indicate  that  the  curvature  modification  we 
have  incorporated  results  in  changes  in  the  desired  direction. 

The  computed  contour  of  the  curvature-dependent  c  ,,  is  seen 
in  Figure  18.  It  is  clear  that  the  constant  value  of  0.09  for 
c u  is  valid  only  in  a  limited  region  of  the  flowfield  (and 
outside  the  recirculation  region).  To  the  right  of  the  c u  =  0.09 
contour,  there  are  local  regions  of  higher  values  of  c  y  (up  to 
0.18).  Inside  the  recirculation  region,  however,  the  constancy 
of  cy  is  seen  to  break  down  dramatically.  We  note  that  the 
correction  leads  to  a  significant  reduction  in  c u  in  the  vicinity 
of  the  separated  streamline  (i.e.,  in  the  curved  shear  layer 
bordering  the  recirculation  region).  The  value  of  cy  =  0.025  was 
(as  noted  in  Paragraph  2.3)  an  arbitrarily  chosen  lower  limit 
in  the  computations,  without  which  c u  would  become  zero  or  even 
negative  within  the  shaded  region  (where  the  concept  of  local 
equilibrium  of  turbulence,  on  which  the  curvature  correction  is 
based,  is  unlikely  to  remain  valid).  Our  computed  c  u  contour  is 
quite  similar  to  that  obtained  for  the  ax  isymmetr ic ,  unconfined 
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Figure  18.  Distribution  of  Curvature  Corrected  c 


annular  jet  of  Reference  33,  except  for  the  nore  complicated 
nature  of  the  contour  near  the  centerline  which  results  from  the 
presence  of  the  central  jet  in  our  configuration.  Note  that  the 
centerline  location  near  the  forward  stagnation  point  is 
characterized  'nv  the  low  c  y  region.  This  appears  to  confirm  our 
suspicions  that  the  underprediction  of  the  forward  stagnation 
point  is  traceable  to  too  high  a  c y  value  (see  Figure  3)  employed 
there  under  the  constant  c ,,  model. 

The  centerline  stagnation  point  results  with  the  curvature 
correction  are  seen  in  Figure  17.  For  the  two  data  points  (viz., 
CO2  flow  of  4  and  8  kg/hr)  shown  in  the  figure,  the  effect  of 
curvature  correction  is  somewhat  masked  by  the  increase  in  the 
inlet  length-scale  parameter  X  from  0.3333  to  0.5556.  However, 
the  effect  of  increasing  X  (without  curvature  correction)  is  to 
move  monotonically  both  the  forward  and  rear  stagnation  points 
upstream.  This  results  in  much  greater  underprediction  of  the 
forward  stagnation  point.  When  the  curvature  correction  is 
introduced,  the  forward  stagnation  point  is  moved  farther 
downstream  and  closer  to  the  measured  value.  Clearly,  the 
prediction  with  X  =  0.5556  is  better  than  that  with  X  =  0.3333. 

In  the  case  of  the  centerline  rear  stagnation  points,  the  effect 
of  the  increase  in  X  and  the  curvature  correction  is  to  reduce 
the  degree  of  overprediction  and  to  lead  to  better  agreement  with 
measured  values.  This  observation  becomes  clearer  from  an 
inspection  of  Figure  19  which  presents  the  centerline  decay  of 
the  mean  axial  velocity.  For  the  case  of  3  kg/hr  CO2  flow,  we 
see  that  the  mean  axial  velocity  profile  shows  better  agreement 
between  prediction  and  measurement  with  curvature  correction  (for 
X  =  0.5556)  than  without  it.  Also  noticeable  is  the  reduced 
underprediction  of  the  forward  stagnation  point  and  the  reduced 
overprediction  of  the  rear  stagnation  point  due  to  the  introduc¬ 
tion  of  the  curvature  correction.  Of  course,  in  the  reverse-flow 
region,  the  uncorrected  profile  seems  to  be  in  better  agreement 


96 


EXPERIMENT,  16  kg/hr  <rn  ,„of 

EXPERIMENT,  B  kg/hnC02* A,r '2k9/s  (Ref 
EXPERIMENT,  4  kg/hr 
COMPUTATION  w/  CURVATURE  CORRECTION 
COMPUTATION  w/o  CURVATURE  CORRECTION 


X1 


1 


i 

i 

i 


with  measurements  in  the  vicinity  of  z/I>0.6.  This  emphasizes 
the  ad  hoc  nature  of  the  curvature  correction. 

The  other  noteworthy  feature  of  Figure  19  is  the  excellent 
agreement  between  the  measured  and  predicted  profiles  for  the 
CC>2  flow  of  16  kg/hr.  As  seen  in  Figure  17,  this  flow  rate  is 
predicted  to  result  in  the  occurrence  of  centerline  reverse  flow 
(with  both  stagnation  points  being  present)  for  X  =  0.3333  and 
without  curvature  correction.  However,  in  Figure  19  it  is  seen 
that  for  X  =  0.5556  and  with  correction  for  curvature  effects, 
the  16  leg /hr  CO2  flow  just  achieves  "breakthrough"  of  the 
recirculation  region.  Since  the  experiments*-®  implied  a 
"breakthrough"  flow  rate  of  14.7  kg/hr,  the  curvature-corrected 
results  exhibit  better  agreement  with  the  experimental  data  than 
did  the  uncorrected  predictions.  Indeed,  the  present  LDA  results 
(e.g.,  see  Figure  32  of  Reference  47)  show  that  the  centerline 
flow  reversal  is  still  present  at  the  CO2  flow  rate  of  16  kg/hr. 
Thus,  the  earlier  implication  that  the  "breakthrough"  occurred  at 
a  flow  rate  of  14.7  kg/hr  is  incorrect.  In  the  light  of  the 
recent  experimental  data,  it  also  appears  that  the  "standard"  k-e 
model  (i.e.,  without  accounting  for  streamline  curvature  effects) 
predicts  the  APL  flowfield  better  than  the  comparison  with  the 
earlier  measurements^  had  indicated.  This  becomes  clear  in  the 
following  discussion. 

4.5  COMPARISON  OF  PREDICTIONS  WITH  THE  NEWER  EXPERIMENTAL 

RESULTS 

In  the  previous  paragraphs  the  influence  of  the  differencing 
schemes,  inlet  turbulence-length  scales  and  streamline-curvature 
effects  on  the  numerical  predictions  was  discussed.  Before  we 
conclude  this  section,  the  numerical  predictions  are  examined  in 
terms  of  the  recent  APL  data  on  the  velocity1*7  and  CO2 
concentration1*®  fields. 
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4.5.1  Centerline  Variation  of  the  Mean  and  rms  Axial  Velocity 

fields  “ 

The  normalized  profiles  of  the  centerline  mean  and  r-ns 
axial  velocity  components  for  three  different  CQ-?  flow  rates  are 
shown  in  Figures  20  through  22.  The  predictions  were  based  upon 
the  (41  x  34)  computational  grid,  the  "hybrid”  upwind  dif¬ 
ferencing  scheme  and  the  "standard"  k-e  model.  The  value  of  X 
was  0.5556  for  the  case  of  zero  COg  flow  and  0.3333  for  the  other 
two  cases. 

The  experimental  data  for  the  mean  axial  velocity  for  the 
CO2  flow  rates  of  0,  6  and  16  kg/hr  differ  from  the  previous 
results.  1-0  The  locations  of  the  rear  stagnation  point  are  at  a 
z/D  of  1  in  Figure  20,  slightly  greater  than  1  in  Figure  21  and 
slightly  less  than  1  in  Figure  22.  The  measurements  in  Figure  22 
for  16  kg/hr  also  indicate  that  for  0.75  <  z/D  <  0.975  the 
centerline-flow  reversal  is  present.  The  earlier  results , 1° » 
on  the  other  hand,  indicated  that  the  rear  stagnation  point  was 
at  a  z/D  of  0.9  (e.g.,  see  Figure  16)  and  that  the  minimum 
centerline  mean  axial  velocity  at  16  kg/hr  was  a  finitely  large 
positive  value  (e.g.,  see  Figure  19).  Thus,  the  newer 
experimental  results  show  better  agreement  with  the  present 
predictions  and  those  of  Reference  6. 

The  agreement  between  the  measurement  and  prediction  for 
the  mean  velocity  is  generally  good  upstream  of  the  rear  stagna¬ 
tion  point,  except  for  the  underprediction  in  Figures  20  and  21 
between  the  location  of  the  peak  negative  velocity  and  the  rear 
stagnation  point.  Because  of  the  very  small  region  of  flow 
reversal  in  Figure  22,  there  is  much  less  disagreement  at  the 
highest  COg  flow  rate.  Indeed,  a  comparison  of  Figures  19  and  22 
reveals  that  the  predictions  with  and  without  curvature  effects 
at  16  kg/hr  do  not  differ  significantly.  This  should  not  be 
surprising  since  the  central  jet  essentially  "breaks  through"  the 
recirculation  region  and  consequently  does  not  suffer  the  large 
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streamline  curvature  effects  associated  with  the  recirculating 
flow.  Finally,  the  measured  recovery  of  the  mean  axial  velocity 
downstream  of  the  rear  stagnation  point  is  greater  than  that 
given  by  the  prediction  for  all  the  three  CO 2  flow  rates.  This 
feature  seems  to  be  characteristic  of  the  "standard"  k-e  model, 
presumably  reflecting  that  the  isotropy  assumption  in  the  model 
may  be  invalid  and  that  the  curvature  effects  may  be  significant. 
This  difference  between  the  measured  and  predicted  recovery  rate 
is  also  consistent  with  the  difference  noted  between  the  measured 
and  predicted  rns  axial  velocity  component  in  Figures  20  and  21. 
The  experimental  results  show  a  pet'  ing  and  a  subsequent  ""ore 
rapid  decay  of  the  turbulence  intensity,  while  the  predicted 
results  indicate  a  continuously  decreasing  trend  at  a  slower 
rate  . 

Figure  21  shows  that  the  measured  forward  stagnation  point 
occurs  at  z/D  =  0.28.  Although  the  calculation  still  underpre¬ 
dicts  this  location,  the  recent  measurement  shows  better 
agreement  with  the  present  predictions  and  those  of  Reference  6. 
For  example,  according  to  the  earlier  measurement  1°  (see  Figure 
16),  the  forward  stagnation  point  at  6  kg/hr  occurred  at  a 
z/D  ~  0.4  which  indicated  considerable  underprediction  by  the 
calculations.  The  much  closer  agreement  between  the  predictions 
and  the  recent  measurement  may  be  attributed,  at  least  in  part, 
to  the  crucial  difference  in  the  central-jet  exit  configuration 
in  the  two  measurements.  The  recent  experiments  employed  a 
we  1 1 -des  igned  nozzle  with  a  contraction  ratio  of  2.56  (e.g.,  see 
Figure  16  in  Part  I  of  this  report),  thereby  ensuring  a  nearly 
uniform  exit-velocity  profile  and  thus  conforming  more  closely  to 
the  assumed  uniform  profile  in  the  calculations.  The  earlier 
experiments^  ,  however,  involved  a  straight  tube  15  liameters  in 
length  upstream  of  the  exit  plane,  thereby  resulting  in  a 
nonuniform  exit-velocity  profile.  Clearly,  the  latter 
configuration  leads  to  a  larger  centerline  axial  velocity  t^an 
the  former  (under  identical  mass  flow  rates)  and  consequently  to 
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a  larger  value  for  the  location  of  the  forward  stagnation  point. 
While  this  difference  in  the  central-jet  geometry  accounts  for 
the  appreciable  decrease  in  the  location  of  the  forward 
stagnation  point,  it  is  unlikely  to  have  contributed  to  the 
observed  increase  in  the  location  of  the  rear  stagnation  point. 

The  quantitative  agreement  between  the  measurement  and 
prediction  for  the  rms  axial  velocity  in  Figures  20  through  22  is 
generally  poorer  than  for  the  mean  velocity.  For  all  three  CO2 
flow  rates,  the  predictions  show  the  observed  trend  (except  for 
the  peaking  near  the  rear  stagnation  point  in  Figures  20  and  21). 
For  example,  in  Figure  21  we  see  that  boch  measurement  and 
prediction  show  the  sharp  peak  in  the  rms  component  in  the 
vicinity  of  the  forward  stagnation  point.  Likewise,  in  Fiqure  22 
there  is  reasonable  agreement  for  the  location  of  the  sharp  peak 
(indicating,  perhaps,  the  transition  to  turbulence)  between  the 
prediction  and  measurement.  As  noted  earlier,  the  assumed 
isotropic  turbulence  and  lack  of  accounting  for  curvature  effects 
in  the  turbulence  model,  among  other  things,  suggest  that  the 
agreement  between  the  prediction  and  measurement  of  the 
turbulence  structure  in  the  centerbody  flowfield  is  unlikely  to 
be  as  good  as  that  noted  for  the  mean  velocity  field. 

4.5.2  Centerline  Variation  of  COg  Concentration 

Figures  23  and  24  show  the  comparison  of  the  predicted 
centerline  mole  fraction  of  COg  with  the  APL  (intrusive)  gas 
sampling  measurements . 43  The  agreement  bewteen  the  prediction 
and  measurement  for  the  overall  trend  in  both  COg  flow  rates  is 
good.  The  quantitative  agreement  for  16  kg/hr  is  better  than 
that  for  6  kg/hr.  This  behavior  appears  consistent  with  that 
noted  for  the  mean  axial  velocity  in  Figures  21  and  22. 

We  note  that  the  earlier  anticipated  trends^  of  the  cen¬ 
terline  COg  concentrations  showing  a  rapid  decay  first  and  an 
equally  rapid  approach  to  uniform  values  subsequently  in  the 
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annular-flow  dominant  condition  have  been  confirmed  by  the 
present  predictions  and  measurements . 48  Furthermore,  as  can  be 
seen  by  comparing  Figures  21  and  23,  and  22  and  24,  the  points  of 
intersection  obtained  by  extrapolating  the  portions  of  the 
profiles  that  denote  the  rapid  decay  and  the  approach  to 
unifornization  in  Figures  23  and  24  fall  very  close  to  the 
forward  stagnation  points  of  Figures  21  and  22  respectively.  In 
view  of  the  rationale  presented  in  Reference  5  for  expecting  the 
turning  of  the  centerline  concentration  profile  to  occur  in  the 
vicinity  of  the  forward  stagnation  point,  the  internal 
consistency  demonstrated  by  both  the  predictions  and  measurements 
is  gratifying. 

4.5.3  Comparison  of  the  Measured  and  Predicted  Streamlines 

It  is  instructive  to  examine  the  time  averaged  contours4^ 
of  the  normalized  stream  function  shown  in  Figure  25.  The 
normalization  is  with  respect  to  the  inlet  annular  air  mass  flow 
of  2  kg/s.  The  contour  labelled  11  denotes  the  zero  stream- 
function  contour  which  separates  from  the  centerbody  surface  and 
meets  the  centerline  at  the  rear  stagnation  point.  The  contours 
(outside  the  separated  streamline)  labelled  1  through  10  denote 
the  (nonrecirculating)  annular  stream.  The  closed  contours 
labelled  12  through  18,  confined  between  the  centerline,  the 
centerbody  face  and  the  zero  stream-function  contour,  represent 
the  time-averaged  recirculating  vortex. 

The  upper  part  of  Figure  25  corresponds  to  the  measured 
contours,  reproduced  from  Figure  30  of  Part  I  of  this  report. 
Unlike  the  smooth  contours  of  Figure  30,  the  contours  in  Figure 
25  were  obtained  by  AFWAL/POSF  through  computer  graphics. 49  The 
lower  part  of  Figure  25  corresponds  to  the  predicted  contours 
obtained  from  the  TEACH-Code  computations  performed  on  the 
MODCOMP  computer  system  at  ^FWAL/POSF  by  Mr.  J.  S.  Stutrud.^9 
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The  comparison  of  the  measured  and  predicted  contours 
reveal  several  interesting  features.  The  axial  extent  of  the 
recirculation  region  is  slightly  underpredicted.  This  conforms 
to  the  earlier  experience  of  several  investigators.  There  is  a 
progressive  degree  of  underprediction  of  the  reverse-mass  flow 
contours.  The  strength  of  the  vortex  expressed  in  terms  of  the 
peak  negative  mass  flow  is  measured  to  be  6.5%  of  the  inlet  mass 
flow  (the  contour  labelled  13  is  taken  to  represent  the  vortex 
center  which  corresponds  to  the  negative  maximum).  The  predicted 
value  is  about  5.1%  of  the  inlet  mass  flow.  The  FREP-Code 
calculations  (which  employed  a  constant  eddy-viscosity  model)  in 
Reference  5  predicted  the  strength  of  the  vortex  to  be  5.37%. 

This  indicates  that  the  "standard"  k-e  model  tends  to 
underpredict  the  strength  of  the  vortex.  Since  Figure  19  shows 
that  the  prediction  accounting  for  streamline  curvature  effects 
results  in  a  higher  peak  negative  axial  velocity  (and 
consequently  a  higher  vortex  strength),  it  would  seem  that  the 
interior  of  the  recirculating  region  is  significantly  influenced 
by  curvature  effects. 

The  above  conclusion  is  also  borne  out  by  an  examination 
of  the  location  of  the  vortex  center.  The  measurement  in  Figure 
25  indicates  that  the  axial  and  radial  coordinates  (normalized  by 
the  centerbody  diameter  D)  of  the  vortex  center  are  respectively 
9.43  and  0.35.  The  predicted  results  from  Figure  25  are  0.26  and 
0.35.  The  predicted  values  reported  in  Reference  5  are  0.3  and 
0.35.  The  very  good  agreement  seen  in  Figure  25  between  the 
predicted  and  measured  radial  coordinate  is  surprising, 
especially  when  there  is  significant  underpred ic t ion  of  the  axial 
coordinate.  Also,  the  identical  value  predicted  by  the  constant 
eddy-viscosity  model^  and  by  the  more  realistic  k-e  turbulence 
model  here  implies  that  the  details  of  turbulence  model  do  not 
have  much  effect  on  how  far  the  vortex  center  is  radially 
displaced  from  the  centerline.  Indeed,  as  noted  in  Reference  5, 
the  radial  coordinate  in  terms  of  the  duct  diameter  is  0.19  and 
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this  is  quite  to  the  value  of  0.13  reported,  hy  '<o  and 

Chan  9  ■*  for  their  ancon f  i  neh  annular  jet.  In  view  of  their 
observation  fat  the  racial  position  of  the  vortex  center  is 
essentially  m-ieoen  d~ "  t  of  the  noment  jn  flux  of  the  annular  ^et 
( and  thus  the  pressure  available  for  entrainment  behind  the 
centerbody  face),  it  appears  that  the  radial  coordinate  is 
determined  more  by  the  geometry  than  hy  the  flowfield  details. 

The  disagreement  noted  for  the  axial  coordinate  (between  the 
measurement  and  orediction  in  Figure  25  as  well  as  between  the 
previous^  and  present  predictions),  however,  indicates  that  the 
axial  coordinate  does  depend  on  the  flowfield  structure.  A 
possible  line  of  speculation  in  this  regard  is  suggested  by  the 
results  of  Ko  and  Chan'*'*  which  show  that  the  axial  coordinate 
decreases  with  increasing  nond imens iona 1  momentum  flux  of  the 
annular  jet  and  consequent  increasing  entrainment  behind  the 
centerbody  face.  The  present  discussion  is  not  concerned  with 
the  variation  in  the  annular  momentum  flux.  Nevertheless,  we 
speculate  here  that  the  predictions  overestimate  the  entrainment 
by  the  annular  jet  as  compared  to  the  measurement.  Likewise,  the 
present  k-e  model  overpredicts  the  entrainment  when  compared  to 
the  eddy-viscosity  model. 5  Extending  this  line  of  speculation 
further,  it  would  seem  that  the  streamline  curvature  correction 
to  the  k-e  model  tends  to  diminish  the  entrainment  and  thereby 
move  the  axial  coordinate  of  the  vortex  center  further 
downstream.  This  conjecture  can  be  verified  when  the  entrainment 
characteristics  of  the  flowfield  are  examined  through  measurement 
and  prediction. 

Before  concluding  this  discussion,  it  is  worth  mentioning 
that  there  are  some  differences  between  the  MODCOMP  predictions*9 
seen  in  Figure  25  and  the  predictions  seen  in  Figures  20  through 
22.  For  example,  the  calculations'*9  employed  a  (41  x  39) 
computational  grid,  a  computational  domain  of  <50  cm  in  the  axial 
direction,  power-law  differencing  scheme  and  a  value  of  0.5556 
for  Accordingly,  Figure  20  for  the  case  of  zero  Cb2  *low 
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shows  a  larger  value  than  Figure  25  for  both  the  axial  coordinate 
of  the  vortex  center  (which  is  close  to  the  axial  location  of  the 
centerline  peak  negative  Tiean  axial  velocity)  and  the  rear 
stagnation  point.  However,  the  underprediction  (with  respect  to 
the  measurement  in  Figure  25)  of  the  axial  coordinate  remains  and 
our  speculative  arguments  still  hold. 
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SECTION  5 

CONCLUSIONS  AND  RECOMMENDATIONS 
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In  this  section  the  major  conclusions  from  the  numerical 
flowfield  modeling  of  the  centerbody  configuration  are  presented, 
followed  by  recommendations  for  further  research. 

5.1  CONCLUSIONS 

a.  The  numerical  modeling  has  successfully  demonstrated  the 
capability  of  the  TEACH-T  computer  program  to  provide  physically 
correct  predictions  of  the  isothermal  turbulent  recirculating 
flowfields  of  the  centerbody  combustor  configuration. 

b.  The  performance  of  the  "standard”  k— ;  turbulence  model 
in  the  teach  code  has  been  free  from  the  numerical  convergence 
difficulties  with  the  FREP  code  in  earlier  modeling  activities. 
Thus,  the  present  modeling  avoids  any  arbitrariness  inherent  in 
the  assumption  of  a  constant  eddy  viscosity  model. 

c.  The  numerical  calculations  of  the  flowfields  in  both  the 
APL  and  UCI  configurations  predict  flowfield  features  that 
conform  to  the  experimental  observations. 

d.  The  present  predictions  of  the  APL  configuration  are  in 
good  agreement  with  the  numerical  results  of  this  flowfield 
reported  in  the  literature, 

e.  In  the  annular-flow  dominant  regime  where  two  stagnation 
points  occur  on  the  centerline,  the  present  predictions  show  much 
better  agreement  with  the  measured  results  of  stagnation  points 
than  do  the  FREP  predictions.  As  in  the  case  of  the  predicted 
data  available  in  the  literature,  the  underpred ict ion  of  the 
forward  stagnation  point  and  the  overprediction  of  the  rear 
stagnation  are  present  to  a  small  extent.  This  appears  to  be  due 
to  the  "standard"  k-e  turbulence  model  which  does  not  include  the 
effects  of  streamline  curvature  and  the  "hybrid"  upwind 
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differencing  scheme  which  retains  diffussive  effects  only  for 
cell  Peclet  numbers  less  than  or  equal  to  2. 

f.  In  the  central-jet  dominant  regime  (when  the  annular  air 
flow  is  very  small  or  when  it  is  completely  absent),  the 
centerline  decay  of  the  mean  axial  velocity  resembles  that  of  a 
free  jet.  Both  the  normalized  velocity  data  and  CO 2  mass 
fraction  data  reduce  to  a  single  curve  which  has  the 
characteristic  slope  of  -1,  when  displayed  on  a  log-log  plot. 

The  distinction  between  small  annular  flow  and  zero  annular  flow 
in  the  FREP  calculations  is  not  observed  in  the  present  work. 

g.  The  present  study  shows  that  the  radial  distributions  of 
the  mean  axial  velocity  field  exhibit  self-similarity  when 
Abramovich-type  normalizations  of  the  axial  velocity  and  radial 
distance  are  employed.  This  behavior  is  noted  for  both  the  UCI 
and  APL  configurations,  for  different  annular  and  central  jet 
flow  rates  (when  either  the  annular  jet  is  dominant,  or  the 
central  jet  is  dominant,  or  neither  jet  is  dominant),  and  for  a 
range  of  axial  locations  (within  and  without  the  recirculating 
region).  This  similarity  is  also  observed  for  the  measured 
velocity  data  in  the  APL  configuration,  the  unconfined  annular 
jet  configuration  (Ko  and  Chan44),  and  the  unconfined  flow 
downstream  of  a  disk  (Durao  and  Whitelaw4^).  Although  the 
similarity  may  remain  incomplete  in  certain  cases  due  to  other 
factors  (such  as  the  nonvanishing  mean  radial  velocity  component 
at  the  inlet  for  the  flow  past  a  disk),  it  does  appear  that  the 
measured  and  predicted  results  exhibit  the  tendency  toward 
similarity.  We  believe  that  this  flowfield  similarity  may  well 
be  a  necessary  condition  for  the  correctness  of  the  measured  and 
predicted  data. 

h.  With  the  nonuniform  grid  spacing  adopted  for  the  41 
axial  x  34  radial  grid-point  distribution  and  the  exit  plane  of 
the  computational  domain  located  at  least  more  than  twice  the 
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centerbody  diameter  from  the  inlet  plane,  the  present  results 
generally  appear  to  be  grid-independent.  Sensitivity  tests 
involving  grid  spacings  that  are  50%  larger  and  exit  plane 
location  that  is  50%  farther  have  resulted  in  centerline 
stagnation  point  predictions  that  vary  by  less  than  10%. 

i.  Sensitivity  tests  for  the  influence  on  the  predicted 
results  of  boundary-layer  thickness  in  the  inlets,  turbulence 
intensity  in  the  inlets,  and  the  Schmidt  number  for  the 
dissipation  equation  do  not  appear  to  result  in  significant 
variations  in  the  parametric  range  tested. 

j.  Sensitivity  tests  for  the  influence  of  the  inlet 
turbulence  length  scales  and  the  k- e  turbulence  model  constant 
cu  show  that  the  predicted  results  of  the  centerline  stagnation 
points  are  greatly  affected  by  the  assumed  values  of  the  inlet 
length  scales  and  cu. 

k.  Modifications  of  the  TEACH  code  implemented  in  the 

a 

present  program  by  the  replacement  of  the  "hybrid”  upwind 
differencing  with  the  power-law  differencing  scheme  result  in 
improved  numerical  predictions,  in  the  limited  parametric  testing 
completed . 

l.  Ad  hoc  modifications  to  the  "standard"  k- t  model  for 
incorporating  the  effects  of  streamline  curvature  through  the 
introduction  of  a  curvature-corrected  and  hence  nonconstant 
cu  appear  to  result  in  predictions  that  show  better  agreement 
with  the  experimental  data  than  the  predictions  with  the 
"standard"  k-e  model. 

• 
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5.2  RECOMMENDATIONS  FOR  FURTHER  ACTIVITY 


Based  on  the  foregoing  conclusions,  we  offer  the  following 
recommendations  for  future  research  in  numerical  modeling. 

a.  For  isothermal  flowfield  modeling  with  the  ""EACH  code, 
additional  improvements  in  the  differencing  schemes  (e.g.,  the 
skew-upwind  differencing  scheme  and  the  quadratic,  upstream- 
weighted  differencing  scheme)  are  possible,  and  this  aspect 
deserves  study. 

b.  In  view  of  the  ad  hoc  nature  of  the  streamline  curvature 
correction  to  the  "standard"  k-  e  model  implemented  in  the  present 
program,  a  more  systematic  approach  and  rigorous  formulation  for 
including  the  effects  of  streamline  curvature  must  be 
investigated . 

c.  Isothermal  flowfield  calculations  of  the  APL  and  UCI 
configurations  must  be  carried  out  with  the  viewpoint  of 
establishing  scaling  criteria. 

d.  Isothermal  flowfield  calculations  of  the  APL 
configuration  must  be  performed  to  establish  the  influence  of 
annular  flow  rates  and  blockage  ratios  on  the  location  of  the 
vortex  centers  and  the  strength  of  the  recirculation  region. 

e.  Modification  of  the  presently  available  version  of  the 
TEACH  code  or  the  use  of  any  refined  version,  if  available,  must 
be  considered  for  the  modeling  of  reacting  flowfields. 

f.  When  the  reacting  flowfield  predictions  are  available, 
the  implication  of  the  vortex  center  on  flame  stabilization  must 
be  studied. 

g.  Time-dependent  flowfield  calculations  must  be  performed 
to  investigate  the  temporal  characteristics  of  the  centerbody 
configuration . 
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